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ABSTRACT 


A review is given of current computational methods for analyzing flows in 
turbomachinery and other related internal propulsion components. The methods 
are divided primarily into two classes, inviscid and viscous. The inviscid 
methods deal specifically with turbomachinery applications. Viscous methods, 
on the other hand, due to the state-of-the-art, deal with generalized duct 
flows as well as flows in turbomachinery passages. Inviscid methods are 
categorized into the potential, stream function, and Euler approaches. 

Viscous methods are treated in terms of parabolic, partially parabolic, and 
elliptic procedures. Various grids used in association with these proced 
are also discussed. 


The subject of internal flows is a very broad and complex one, encompassi' 
wide variety of geometries and flow situations. There are many examples _ 
machinery in which internal flows play an important part. One such example is 
the modern turbofan engine, such as that depicted in Figure 1. Here the 
understanding of internal flows is very important in predicting the perfor- 
mance of many key components. These include the inlets and exhaust nozzles at 
the extremeties of the engine, the rotating and stationary turbomachinery 
blade rows in both the compressor and turbine sections of the engine, the 
interconnecting ducts, and the combustor portion of the engine. Another 
significant type of machinery in which internal flows are extremely important 
is the internal combustion engine. Here the flows to the various inlet and 
exhaust ducts, ports, and cylinders are extremely complex and basically 
unsteady. These complex flows must be understood before the performance of 
such engines can be predicted. Another situation where complex internal flows 
play a large role in the performance of machinery is that of nuclear reactors. 
Here the flows through a maze of tubes and passages interact to influence the 
performance and reliability of the entire system. In addition, the steam 
generator portion of the nuclear reactor has all complexity of the gas turbine 
engine described previously, as well as the complexities brought on by a two 
phase flow situation. 

In this paper, in order to make the subject more manageable, we have chosen to 

treat in detail a subset of the total class of internal flows. We will be 

speaking specifically about flows through turbomachinery blade rows of all 
types, as well as viscous flows through ducts of various geometries. 

The paper will be divided primarily into descriptions of inviscid flow methods 

and then viscous flow methods. Inviscid flow methods will be described in the 
context of turbomachinery blade row applications because a large number of 
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analyses have been developed for these situations. For viscous flow methods 
the state-of-the-art is less well advanced. Here the treatment will be 
expanded beyond turbomachinery to also encompass duct flows. The analyses 
discussed in the paper will also be limited to steady flows. Although some do 
have the capability to predict unsteady flow phenomena, they have been 
developed primarily as predictors for steady flows. They generally are not 
used in their present form in order to solve unsteady flow problems in turbo- 
machinery. 

In describing the analyses, a variety of different geometries will be discus- 
sed. We will be describing analyses which apply first of all to axial blade 
geometries in both compressor and turoine blade rows. Some of these analyses 
are also designed to predict mixed flow situations, which are encountered in 
both centrifugal compressors and radial inflow turbines. There are also some 
analyses designed for pure radial situations, such as those encountered in 
radial diffusers as well as the inlet pre-swirl vanes in radial turbine situa- 
tions. Some of the above analyses also appl> to geometries which have no 
blades, such as the interconnecting ducts and bends in turbomachinery. 

Another differentiation is that between rotating and stationary blade rows, 
with most of the analyses applying to both. Finally, some types of analyses 
apply only in certain flow ranges. Most will handle the subsonic flow regime, 
while emphasis in the recent past has been upon analyses for the transonic and 
supersonic flow regimes. 

A wide variety of flow characteristics exist in the various types of turboma- 
chinery just described. The objective of the analyses to be discussed later 
will be to predict as many of these flow features as is possible. Figure 2 
illustrates for a turbomachinery blade row situation many of the important 
flow features which we will ideally be trying to compute with the methods 
described in this paper. First of all, large axial, radial and centrifugal 
pressure gradients exist within the flow passages due to the turning of the 
fluid within the blade rows. Secondly, this turning of the fluid within the 
passages redistributes the incoming vorticity field and generates cross 
flows. At higher velocities strong shocks exist within the blade passages. 
These can be complex and interacting, and they in turn can generate their own 
additional vorticity fields. These shocks, of course, interact with the blade 
surfaces and endwall boundary layers, often causing separation and additional 
blade loss. 

There are also a wide variety of viscous flow phenomena existing in the blade 
passages. Primarily there are the boundary layers which exist on all blade 
and endwall surfaces as well as on the surfaces of midspan dampers which may 
exist within the blade passages. Such boundary layers can have laminar, 
transitional, ard turbulent regions. When pressure gradients are strong, of 
course, separation can also occur. Some separations may experience reattach- 
ment. Finally, there are the wakes which exist downstream of all blade rows. 

In addition to these common flow features, there are other more complex 
viscous flows which serve to redistribute the internal vorticity in the blade 
passages. The first of these, commonly encountered in turbines, is the horse- 
shoe vortex which is generated at the junction of the blunt leading edge and 
the endwall. Such vortices curl up, flow across the passage, ana are shed 
downstream off the endwall or the surface of the adjacent blade row. Another 
complex flow phenomenon is the tip clearance flow in which a vortical flow is 
inducea by the leakage of fluia across the unshrouded tip of a rotating or a 
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stationary blade. Finally, the fact that half the blade rows operate in a 
rotating reference frame introduces the effects of centrifugal and Coriolis 
forces on both the mean flow and turbulence in these flow passages. 

The resultant flow picture is an extremely complex one, particularly in the 
multi-stage environment that exists in modern turbomachinery. No single analy- 
sis can hope to model all of these flow phenomenon at the same time. In the 
analyses to be described here, a number of different approaches are used to 
divide this overall problem into one of manageable size. We have chosen here 
to separate these models into two major groups, inviscid analyses and viscous 
analyses. Within these classifications there are many other types of modeling 
assumptions which are made. First, the assumption is always made to limit 
these analyses to a single blade row, either stationary or rotating (except in 
the case of hub-to-shroud stream function analyses, some of which have multi- 
stage capability). Next, within a given blade row, for a steady state type of 
solution, the assumption is made that the flow in all blade passages is iden- 
tical; thus only a single blade passage has to be analyzed. Within a blade 
passage the next decision to be made is whether a two-dimensional or three- 
dimensional flow solution will be obtained. Obviously in a three-dimensional 
flow solution the entire passage is considered; however, many times assump- 
tions are made to eliminate flows over the blade tips or the consideration of 
dampers which may exist in the blade passage. Two-dimensional flow solutions 
are typically done on one of two types of flow surfaces. The first of these 
is a meridional plane solution, treating equations that model flow on an 
average or mean flow surface between the suction and pressure sides of the 
blade (Figure 3). This surface is generally formed from the projection of the 
mean camber line of the blade onto the mid passage, with corrections made at 
the leading and trailing edges for incidence and deviation. The second type 
of surface on which 2D analyses can be made is the blade-to-blade surface. 

This surface is generally formed by rotating around the blade row, in the 
circumferential direction, one of the calculated streamlines in the meridional 
plane. Such surfaces, therefore, are axisymmetric and usually have radius 
change from inlet to outlet. Quasi-3-dimensiona1 effects can also be consi- 
dered on such a surface by giving it a thickness which varies in the meridio- 
nal flow direction from inlet to outlet. 

Beyond the decision concerning the dimensionality of the solution, a number of 
decisions can be made in the process of modeling the full flow equations down 
to a reasonable subset to be solved in the particular analysis being performed. 
Tradeoffs have to be made between modeling sufficient flow physics to capture 
the important features of the flow, and reducing the equations and boundary 
conditions such that exorbitant computer time is not consumed in obtaining a 
solution. Some of the assumptions and approaches which have been used in the 
past will be summarized in the historical section which follows. 


Historical Perspective 

A number of excellent articles have appeared in the literature in the last 
several years reviewing the different analysis methods and theories which have 
been used to describe the fluid flow through turbomachinery [1, 2, 3, 4, 5, 6, 
8]. These do an excellent job of reviewing the early work in the field, as 
well as some of the more recent approaches which have been developed. 

In the 50' s and 60' s, singularity methods were often used to compute the two- 
dimensional incompressible potential flows through cascade blade rows. Such 
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analyses were primarily applied on blade-to-blade types of flow surfaces. 

These analyses have been reviewed by Gostelow [1], and Perkins and Horlock 
[3]. These methods will not be discussed in further detail here. 

Also in the 50's and 60's classical secondary flow theory was developed to 
predict three dimensional incompressible rotational or vortical flow in cas- 
cades. Such methods have been reviewed by Horlock and Lakshminarayana [2]. 
These methods likewise will not be discussed in this paper. 

In the 60's and early 70's finite difference approaches began to gain promi- 
nence, and were used to calculate two-dimensional, subsonic, inviscid flows in 
turbomachinery. Both streamline curvature and stream function approaches were 
used, and were applied on both meridional and blade-to-blade flow surfaces. 
These approaches are reviewed by Gostelow [1], Perkins and Horlock [3], and 
Japikse [4]. Stream function methods will be discussed in the inviscid flow 
sections to fol low. 

During the 1970's time-marching solutions of the Euler equations began to be 
used to solve both two-dimensional and three-dimensional transonic flows in 
blade rows. In the 2D cases such methods were applied on both meridional and 
blade-to-blade surfaces. Some of these methods are reviewed by Japikse [4], 
and Habashi [5]. They will be discussed and updated in this paper. 

Also during the 1970' s numerical solutions to the full potential equations for 
two dimensional transonic flow in turbomachinery began to appear. A great 
deal of work has been done to extend such methods up to the present time. 

Early full potential solutions were reviewed by Habashi [5], and the later 
approaches will be discussed and updated here. 

During the 70's many of the above methods, as well as some early viscous 
approaches, began to be applied to flows in centrifugal impellers. These 
methods are reviewed in the paper by Adler [6]. These methods are also 
discussed and updated here. Finally the state-of-the-art in the computation 
of a wide variety of turbulent flows was addressed at the recent 
AFOSR-HTTM-Stanf ord Conference on Complex Turbulent Flows [7]. The relevant 
findings from this conference will also be discussed. 

The review to be given here begins with the consideration of computational 
meshes, since to a large degree the success of analysis approaches hinges on 
the nature of the solution grid which is used. This discussion will include 
consideration of desirable grid properties, the various mesh topologies in 
use, and methods of grid generation. 

Inviscid flow analysis methods are considered next, the Euler equations are 
introduced, and the difficulties in solving the primitive variable form are 
discussed. The stream function formulation for two dimensional flows and the 
scalar potential approximation are both presented, and the advantages and 
limitations of each are described. The various methods currently in use for 
turbomachinery flow analysis are then reviewed. 

The next section will consider viscous methods. The time-averaged Navier- 
Stokes equations are introduced and the many difficulties associated with 
their solution are discussed. Parabolic and partially parabolic approxima- 
tions are presented and the advantages and limitations of each are described. 
The various methods in use for both turbomachinery blade row analyses and for 
generalized duct flow analyses are then reviewed. 
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Finally a short discussion on the status of turbulence modelling is given. 
The paper will conclude with comments on future directions. 


COMPUTATIONAL GRIDS 

The general approach to numerical solution involves discretizing the fluid 
equations on a network of points, or grid, throughout the physical domain. 

The accuracy of the resultant solution can depend to a great degree on the 
properties of this grid. 

Cartesian grids, when used, have generally been used in the physical domain. 
However, they often create problems. At boundaries not aligned with one of 
the particular mesh directions, either partial mesh cells must be used or 
external dumny grid points must be introduced. This results in a cumbersome 
treatment of the boundary conditions and inaccurate results in some cases. A 
second problem is the lack of resolution of the grid in regions of high flow 
gradients. Adequate resolution of the flow in regions adjacent to highly 
curved boundaries often requires local subgrids to be used. The placement of 
subgrids within the larger grid structure necessitates increased coding, and 
the resultant logic problems that come with the interaction between the sub- 
grid and the major grid. A final problem is that boundary layers and wakes, 
which are not generally aligned with the major directions, are often difficult 
to resolve. 

Recent efforts have been placed on devising methods of generating grids 
possessing a number of desirable features which provide a much more accurate 
solution in the important regions of the flow. Some of the characteristics of 
such grids are the following. 

1. At solid surfaces all grid lines shoulo be either parallel or normal to the 
boundary. This simplifies the boundary condition formulation and thus 
improves the accuracy of the solution. It can also simplify the modelling 
of the viscous terms in viscous flow calculations. 

2. Away from solid surfaces the mesh should not be highly sheared, as this 
tends to degrade the accuracy of the solution, especially near shocks and 
other high flow gradient regions. 

3. The grid should provide the ability to cluster lines to a desired density 
in critical regions such as near solid surfaces, at shocks, near leading 
and trailing edges, and in other regions of high surface curvature. 

4. In steady state turbomachinery blade row calculations the assumption is 
usually made that the flow is identical in neighboring passages. 

Therefore mesh periodicity is a requirement to simplify the application of 
boundary conditions at these locations. 

5. At upstream and downstream boundaries, adequate but not excessive mesh 
density is required in order to resolve the flow at these locations and 
properly provide for the application of boundary conditions. 

6. For numerical stability, some algorithms require mesh aspect ratios which 
don't depart too far from one. 
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It is usually not possible to obtain a grid which incorporates all of the 
above features. In many applications these requiren.ents can be in direct 
conflict. For example, the requirement for local clustering of the grid may 
lead to excessive shearing in certain regions. Or, the periodicity requirement 
can be incompatible with the requirement for near-orthogonal ity in staggered 
blade rows. Since no grid can possess all of the above features, obviously 
compromises have to be made which depend on the particular application. 

To satisfy these requirements, several types of grids have been developed for 
use with different geometries. For duct flow calculations, so-called channel 
grids are used (Figure 4). These consist of a set of throughflow lines toge- 
ther with a set of lines or surfaces normal to these. 

For turbomachinery blading three types of grids have been developed. 

1. H-type, which is simply a channel grid which has throughflow grid lines 
aligned with the blades, and the opposing set of grid lines running in the 
circumferential direction (Figures 5 and 6). 

2. 0-type, which has one family of lines forming closed loops about the 
blades, and a second family crossing the first with grid lines radiating 
out from the blade surfaces to the periodic boundaries (Figures 7 ano 8). 

3. C-type, which has one family of grid lines originating at the downstream 
boundary, progressing upstream parallel to the flow, circling the blade 
and continuing downstream again parallel to the flow. The second family 
crosses the first, originating at the blade and wake centerline and 
progressing out to the periodic boundaries (Figures 9 and 10). 

Grid Generation Methods 


Three general approaches have been employed to generate such grids: conformal 
mapping, differential equation techniques, and algebraic techniques. Recent 
updates on a variety of techniques employing these methods are given in 
Reference 9. 

The conformal mapping approach proceeds generally as follows. The physical 
domain is transformed into a simpler one using complex variable techniques. A 
rectangular or near orthogonal grid is then constructed in the transformed 
region by simple algebriac procedures. The transformation of this simple grid 
back into the physical domain then produces the desired curvilinear mesh in 
that domain. This procedure is relatively rapid, but control over the place- 
ment of mesh points in somewhat limited. The method only lends itself to the 
generation of two-dimensional grids, but these in turn can be stacked to give 
a three-dimensicnal grid. General conformal mapping techniques are surveyed 
in Reference 10. 

In the differential equation technique, the curvilinear coordinates in physi- 
cal space are defined to be the solution of a set of coupled Poisson equa- 
tions, one for each of the coordinate directions ( ,yv ,tr , in three 
dimensions). These equations, however, are somewhat difficult to solve in the 
physical space. They are therefore transformed to a quasi-linear elliptic 
system of equations in the rectangular computational space. In this set of 
equations the cartesian physical coordinates (x, y, z) are the unknowns. This 
system is then solved by standard finite difference techniques (relaxation 
methods) to obtain the locations of grid points (x, y, z) in the physical space. 
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This procedure is quite flexible, and allows considerable control over the 
placement of grid points. However, it is more time consuming than other 
approaches. A good description of this technique is contained in Reference 11. 

Finally, in the algebraic techniques, curvilinear coordinates are constructed 
in the physical space through the use of algebraic interpolating functions. 
Opposing bounding surfaces are represented as algebraic functions of some 
parameter, . These opposing surfaces form the first and last members of 
one family of coordinate curves. The second family of coordinate curves is 
formed by connecting points of equal ^ with interpolating functions in the 
common parameter Y! . The remaining members of the first family of coordinate 
curves are formed by connecting points of equal . This method is poten- 
tially the fastest and the most flexible of all tne grid generation techni- 
ques. However, the greater flexibility does require more interaction from the 
user, and the technique may be slightly difficult to master at first. 

Reference 12 describes one of the more advanced algebraic grid generation 
techniques. 


INVISCID METHODS 
Euler Equations 

The ultimate equations we would like to solve in any internal flow situation 
through ducts or turbomachinery are the Navier-Stokes equations which govern 
general viscous fluid flows, however, solving these equations in their full 
form on modern day computers is still quite time consuming. 

A great deal of information can be gained in most flow situations through 
solving a simplified form of the Navier Stokes equations, in which all of the 
viscous terms have been neglected. These are the governing equations for 
inviscid flows known as the Euler equations. These equations in vector form 
are as follows: 



and f and ^ are the static density and pressure, ^nd O' the Cartesian 
velocity components in the x and y directions, ana ^ is the total energy per 
unit volume. The total energy, , is relateo to the internal energy per 
unit mass, ^ , by the following relation 

^ ' I X. 

The internal energy, g , can be related to the pressure and density by the 
ideal gas relation ^ 


( 3 ; 
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where denotes the ratio of specific heats. 

These equations are presented in terms of the conservation variables of mass, 
momentum, and total energy, and these are commonly the terms for which the 
equations are solved. The equations are also presented in terms of two- 
dimensional cartesian components of velocity for the sake of simplicity. Most 
of the cormients made in terms of the equations to be presented in this paper 
can be extended easily to apply to three-dimensional flows, unless otherwise 
noted . 

In most instances the solution desired is that to the steady state Euler equa- 
tions, with the time terms not included. This system of equations is a 
first-order system of partial differential equations. That fact is the source 
of much of the difficulty in obtaining solutions to these equations in their 
present form. The equations change in form depending upon the local flow 
regime. In subsonic flow their character or type is elliptic, whereas in 
supersonic flow regimes they are of hyperbolic type. This leads to major 
difficulties in numerically solving the equations, since the locations of the 
particular flow regimes are not known a priori. In a totally supersonic flow 
field, where the equations are hyperbolic everywhere, some very efficient 
methods exist for their solution. Employing the method of characteristics or 
a simple marching procedure are two common approaches used with supersonic 
flows. In subsonic domains, however, no good general solution method has been 
devised for solving this system of elliptic first order equations. One cormion 
approach that has been devised for solving these equations in either a sub- 
sonic flow field or a mixed subsonic-supersonic flow field is to reintroduce 
the time terms to the equations. The resultant set of first order partial 
differential equations is hyperbolic throughout the flow field. A steady- 
state solution can be obtained by marching these equations from some initial 
guessed flow field through time until an asymptotic steady-state solution is 
achieved. 

Over the past fifteen years, a large number of algorithms have been developed 
for the solution of the time-dependent Euler equations to a steady state. The 
principal disadvantage in approaching a steady state solution in this manner 
lies in the long computational times which are required to wash out of the 
flow field all of the initial transients which are introduced due to the 
initial assumed solution. These initial conditions generally have included 
within them large perturbations away from the final solution. These give rise 
to perturbation waves which move through the solution domain as the solution 
progresses in time. Steady state convergence is reached only when all of 
these perturbations have completely died out. The Euler equations with the 
time terms included do not inherently have any dissipative effect. Therefore 
the discretization of the equations is the only source of damping introduced, 
and this source is rather small. Generally, to reach a steady state requires 
a large number of iterations, and thus a long computational time. 

A second problem inherent in numerically solving the primati ve-variable Euler 
equations is the fact that storage on the computer is quite excessive. In 
these equations there are four primative variables in 2D, five in 3D. Whether 
they be the density, velocity components, and energy, or the more common 
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conservation variables, there are still four or five to be stored at each grid 
point. This is a major disadvantage compared to alternative methods, to be 
described shortly, in whicn only a single parameter must be stored at each 
grid point. 

A third and final source of difficulty in solving the Euler equations is in 
the boundary conditions that are required by this system. Essentially boundary 
conditions are required for all of the primative variables on each of the 
boundaries. Some natural boundary conditions are supplied easily by the 
physics of the problem; however, other auxiliary boundary conditions must be 
obtained, generally through application of the method of characteristics at 
the boundary. Typically this leads to a rather cumbersome treatment, and in 
most methods approximations are made to the full characteristic relationships. 
As an example we generally require some sort of compatibility relation obtained 
from the method of characteristics to get the pressure at solid surfaces. This 
relationship is often replaced by a more approximate relation which can adver- 
sely effect the stability and accuracy of the overall numerical scheme. 

A number of approaches to solving the full set of Euler equations will be sur- 
veyed later in this paper. The conventional approach to circumventing the 
problems related to solving the full Euler equations for steady inviscid flow 
is to define either a stream function or a potential function and solve the 
resultant second order equations which arise from these definitions. 

Stream Function Equation 

The stream function equation is derived by postulating that the mass flow 
components, Pu_ and ^iT , are obtainted from a scalar stream function as 
follows ' 


where is the stream function. Substitution of these definitions into the 
continuity equation (1) shows that it is automatically satisfied. 

If we introduce these terms into the definition of the vorticity 




( 5 ) 



( 6 ) 


we obtain the following 



— 0 > ( 7 ) 


where 6:^ is the vorticity component normal to the plane of solution. 

To work with this equation we must now obtain relations for the density, ^ , 
and the vorticity, , in terms of the stream function, . From the 
perfect gas relation we can obtain the following 


From the 
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where C3L is the speed of sound, S is the entropy, *R. is the gas constant, 
and y- refers to a reference state of zero entropy. The quantities fr* and 
Okr are reference conditions. The entropy S can vary from streamline to 
streamline, but is constant along any given streamline. 

In order to relate the speed of sound to the velocity components, we use the 
equation for conservation of energy as follows: 


^ (u? + u-*-) 

where is speed of sound based on total conditions. 

o i 

Finally the vorticity can be related to total conditions through Crocco's 
relation. 

-¥ Uob I 

dv\. 



(9) 
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where W. denotes a direction normal to the streamlines, and K^is the total 
enthalpy. 


Rearranging we can obtain the following relationship for the vorticity 
terms of density, speed of sound, and entropy. 



J-. A 



in 
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The stream function equation, Eqn. (7), is a second order partial differential 
equation, which allows us to use the extensive experience we have with relaxa- 
tion procedures for the solution of such equations. The stream function 
formulation retains all of the generality contained in the full Euler equa- 
tions, so that it does permit variation in entropy, total pressure and 
temperature throughout the flow field. Likewise the flows which are calcula- 
ted can be rotational and cover the entire flow range from subsonic through 
supersonic flows. However, this formulation does bring with it some inherent 
loss in generality. It is limited to two-dimensional or axisymmetric flows, 
and is also hampered by the fact that the density in the transonic flow regime 
is a double-valued function of the unknown \|; . Solution for in the 
transonic regime is possible, and has recently been obtained by^hafez [ISl for 
an isolated airfoil application. The cascade version of such a development is 
currently underway. A number of different subsonic stream function solutions 
have been obtained for meridional plane and blade-to-blade plane regions in 
turbomachinery applications. These will be surveyed shortly. 


Working in terms of the stream function solves many of the problems cited 
previously for the full Euler equations. The stream function gives a single 
second order partial differential equation, for which robust and well under- 
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stood relaxation solution methods exist. These solutions can therefore be 
obtained in much less computer time *han those for the first oroer Euler 
systems. The computer storage required for this equatiorT is likewise much 
lower, necessitating only the storage of the stream function at each mesh 
point. Finally, the boundary conditions are more natural and are smaller in 
number. For the stream function generally Dirichlet boundary conditions are 
used at all solid surfaces, and either Dirichlet or Neumann boundary condi- 
tions are imposed at all the open surfaces. The resultant problem can be 
solved in at least an order of magnitude less time than the full primative- 
variable Euler equations. 


Full Potential Equation 


Another approach to circumventing the problems inherent in solving the full 
Euler equations is to define a scalar potential function, ^ , such that the 
vector velocity field is everywhere equal to the gradient of that scalar 
potential . 




( 12 ) 


By definition then the curl of that velocity field will be everywhere zero, so 
that such a flow is automatically irrotational throughout the flow field. 

This does place a restriction on the flow, but brings with it the advantage of 
being able to work in terms of a single second-order equation instead of the 
full primative-vari able Euler set. 


If the velocity field is the gradient of a scalar potential then by definition 
the components of velocity u. and iT are related to the potential as follows. 


u.= 





(13) 


where is the scalar potenfial function. Substituting these definitions 
into the continuity equation yields the full potential equation in conservation 
form in two-dimensional Cartesian coordinates. 



(14) 


Refering to Crocco's relationship, equation (11), we see that for an irrota- 
tional two-dimensional flow, if the vorticity is zero, then the flow must be 
everywhere isentropic, and have constant total temperature or enthalpy 
throughout. Through the isentropic relationship, we then relate density to 
the speed of sound as follows. 



Furthermore, knowing that total temperature is conserved throughout yields the 
following relationship between the speed of sound and the velocity components. 
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u*+ o-*-) 


(16) 


Elimination of the density from the conservation form of the full potential 
equation, along with use of the energy relationship given above, leads to the 
following common non-conservation form of the full potential equation. 



Some of the advantages and disadvantages of working with the full potential 
equation are the following. As with the stream function equation, the full 
potential equation permits the user to work with a single second-order partial 
differential equation, for the scalar function <b . Furthermore, it does not 
have the restriction of being limited to two-dimensional flow situations, and 
full 3D flows can be analyzed, as we will see. The two primary disadvantages 
of this approach are that the flow is limited to being irrotational everywhere, 
and it is also isentropic. The isentropic assumption implies that shock waves 
captured in the transonic regime must be limited in Mach number to a value 
less than about 1.3 in order to be accurate. The irrotational ity conoition 
necessitates a uniform incoming flow in two-dimensional flow situations, and 
the condition that equal a constant in the radial direction in three- 
dimens^ional turbomachinery flows. 

The solution of the full potential equation will admit the existence of dis- 
continuities in the flow field. However, these discontinuities are isentropic 
shocks, which do not represent true physical shock waves because they do not 
satisfy the Rankine-Hugoniot jump conditions. However, these shocks will be 
approximately of the proper strength and will exist in the proper position in 
the flow field if the Mach number of the flow approaching the shock is less 
than or equal to 1.3. 

With regard to the irrotational ity condition for full potential flows, the 
question arises as to whether this equation can be applied in the rotating 
reference frame in turbomachinery applications. In this situation, the user 
must recognize that the equations describe an irrotational flow in the 
absolute reference frame, with a solid body rotation imposed on that flow 
field due to the rotation of the wheel. In two-dimensional flows whether the 
blade row is stationary or rotating, if there is no radius change in the flow 
field, then the total pressure and total temperature will be constant through- 
out. If the flow field does have radius change then the rothalpy and entropy 
will be constant along these 2-0 surfaces. The rothalphy, X » is defined in 
terms of the relative velocity components which are expressed as follows 


w, = ^ 

^ >> 2 : 


= 
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where y is the local radius, and , vJ^ , and are the relative 
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velocity components in the axial, radial and tangential directions, 
rothalpy is then defined as 




Wa + \jJ^ ^ \aJ 




-ni 




The 
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where $l is the static enthalpy and is the wheel speed. Similarly in 
three dimensions, in the absolute frame the total pressure and the total 
temperature will be constant throughout. In the relative frame again the 
rothalpy and entropy will be everywhere constant. 

As was the case with the stream function, using the potential formulation 
allows us to work with a single second order partial differential equation and 
obtain all the advantages of the solution methods that exist for that situa- 
tion. The same comments that where made prior with regard to computer storage 
again apply. Finally with regard to boundary conditions, generally for the 
potential function Neumann boundary conditions will exist at all solid sur- 
faces, and either Dirichlet or Neumann boundary conditions will exist at the 
open surface boundaries of the flow field. These are much easier to incorpor- 
ate than the more complicated boundary conditions described earlier for the 
full Euler equations. 

We now describe in detail the analysis that have been devised for the various 
sorts of equations governing inviscid flow situations. Methods for the full 
potential equation will be described first, followed by those for the stream 
function equation, and finally those for the full Euler equations. 


FULL POTENTIAL EQUATION ANALYSES 

A great deal of progress has been made in the last ten years in the develop- 
ment of solutions to the full potential equation for internal turbomachinery 
applications. A major development at the beginning of this period was the 
paper by Murman and Cole [14], in which they demonstrated a way to properly 
account for the domain of dependence in supersonic flow regions by introducing 
a special backward or upwind differencing. This had the effect of stabilizing 
solutions in those regions and permitting the first real transonic flow calcu- 
lations using the potential equation. 

Several major breakthroughs were also made during this period by Jameson. The 
first of these [15] generalized the concepts introduced by Murman and Cole so 
that they applied to the full potential equation in non-conservation form. 
Solutions to the full potential equation in this form now routinely use 
Jameson's rotated difference scheme, which effectively introduces a higher- 
order term which acts like an artificial viscosity in the supersonic region. 
Hafez [16] later introduced the concept of artificial compressibility to 
accomplish somewhat the same objective when the conservation form of the full 
potential is used. In this case the density is evaluated upstream of the 
point at which it is to be applied, which again stabilizes the equation in the 
supersonic zone. Jameson [17] also proposes a methoo for upwinding the 
density when the conservation form of the potential equation is used. 

A number of different authors, applying these numerical techniques, and using 
either finite difference, finite area (or volume), or finite element methods 
to discretize the full potential equation in either conservation or non- 
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conservation form, have devised several excellent, stable, ano accurate 
methods for solving 2D transonic flows on the bl ade-to-bl ade surfaces of 
Ljrbomachinery. Many of these blade-to-blade analyses also incoporate quasi- 
30 effects through radius change and stream channel convergence. We are 
beginning to see a number of full three-dimensional applications of these 
methods, not only to turbomachinery problems, but also to propeller, wing- 
body, and nacelle problems as well. 

Discretizing the Equation 

One of the principal characteristics that differentiates the manner in which 
solutions are obtained to either the conservation or non-conversation form of 
the full potential equation, is the way in which the equation is discretized 
in order to obtain a series of algebraic equations for solution. There are 
primarily three major approaches to this problem: the finite difference 

method, the finite area or volume method, and the finite element approach. 

In the finite difference method the terms in the basic partial differential 
equation, in either untransformed physical space or transformed to a computa- 
tional space, are discretized using standard central differencing and backward 
(or upwind) differencing techniques. The resultant discretized equations, 
written at each of the mesh points, form a linear system of algebraic equa- 
tions which are solved for the unknowns by standard numerical analysis techni- 
ques. This approach was applied in the early work of Dodge [18], Ives [19, 
20], and Rae [21]. More recently this method has been used by Caspar [22, 

23], who applies it to so-called finite or control areas throughout his physi- 
cal domain. 

The next method for discretizing the partial differential equations is termed 
the finite volume method (finite area method in two dimensions). This 
approach is a hybrid between a finite difference method and a finite element 
method, with more of the flavor of the latter. In many algorithms for solving 
the full potential equation, it is transformed and subsequently solved in the 
transformed computational plane. However, this approach has drawbacks, in 
that the transformations can be quite complex, involving a large amount of 
computational labor in evaluating the coefficients of the transformation. The 
desire to have an algorithm which does not depend upon the details of the 
transformation is what motivated the development of finite volume techniques. 
In this method, a local transformation is done on each of the small volumes 
which discretize the flow field, and the transformation metrics are evaluated 
numerically using the Cartesian coordinates of the corners of each of the mesh 
cells. The unknowns in the problem are then represented by some sort of func- 
tional representation on each of these discrete volumes. This method is 
described more fully in the paper by Caughey and Jameson [24]. This method 
has been applied recently in the analyses of Dulikravich [25, 26], Farrell and 
Adamczyk [27], and Fruhauf [28]. 

The third approach is the finite element method. In the finite element method 
the physical space is discretized with a series of triangular or quadrilateral 
shaped elements generated in a completely arbitrary fashion. The generation 
can be readily designed to concentrate elements in regions of high surface 
curvature and where flow gradients are strongest. Once the element grid is 
formed, the potential function is then approximated within each element by 
some sort of linear combination of mesh point values for the function based on 
locally defined shape functions assigned to each mesh point. One of the 
standard weighted residuals methods, usually the Galerkin method, is then used 
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to reduce the partial differential equations to a system of algebraic equa- 
tions which can be solved directly. To apply this method, the continuity 
equation is multiplied by the nodal weighting functions, and is integrated 
over the volumes of the finite elements. The approximation for the potential 
function is then substituted into the continuity equation, and the resulting 
expression represents the error due to the approximation which is being used 
within the elements. The algebraic sum of all such contributions in the 
elements with common nodes are then set equal to zero. This represents the 
algebraic finite element equation for that particular nodal parameter. This 
set of linear equations is solved to yield the potential function throughout 
the flow field. 

The order of accuracy of the functional approximation on each element can be 
enhanced by using more elaborate, higher order elements. Typically linear 
triangular elements are used, but biquadratic and other types of elements can 
be used to give higher order discretization. 

Finite element methods have been used extensively in the blade-to-blade calcu 
lation methods of Laskaris [29, 30], Ecer and Akay [31, 32], and Hirsch and 
Deconinck [33]. The Laskaris work is primarily subsonic, although a three- 
dimensional application is presented in [30]. The Ecer and Hirsch algorithms 
have been extended into the transonic range for turbomachinery applications. 
Eberle [34] has also done a considerable amount of work applying the finite 
element method to transonic potential flow computations. He has applied such 
methods to wings, axisymmetric bodies and nozzles, as well as turbomachinery 
cascades. 

Form of the Potential Equation 


As mentioneo previously, there are two principal forms of the full potential 
equation for which solutions are generally obtained, the conservation form 
(14) and the non-conservation form (17). Each has unique characteristics 
which require that somewhat different methods be applied in their solutions, 
but there doesn't seem to be any distinct advantage in solving one form over 
the other. In fact, the authors cited in this paper are approximately evenly 
divided between the two approaches, with the later approaches favoring solu- 
tion of the conservation form ccupled with artificial compressibility. The 
non-conservation form of the equations has been solved most recently by 
Dulikravich [25], and previously by Ives [19], Rae [21], and Dodge [18]. The 
conservation form is solved by Caspar [22, 23], Farrell [27], Ecer [31, 32], 
and Hirsch [33]. Fruhauf [28] has solutions for both the conservation and 
non-conservation forms of the equation. 

Stabilization in the -Supersonic Zone 


Another characteristic which differentiates among the methods used to solve 
the full potential equation is the technique used to stabilize the equations 
in the supersonic zone for transonic types of calculations. As mentioned 
previously, there are two principal techniques which are used in this regard, 
the artificial viscosity approach of Jameson [15], and the artificial compres 
sibility approach of Hafez [16]. 

The intent of both of these approaches is to modify the differencing in the 
critical flow regions such that the grid points which contribute to the solu- 
tion at a given field point lie within the zone of influence of that point. 
The intent is to obtain finite difference information primarily from grid 
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points which lie within the Mach cone effecting a particular grid point for 

which the equations are being written. Jameson's rotated scheme identifies 
both hyperbolic and elliptic operators in supersonic zones. Upwind differ- 
ences are used for the former and central differences for the latter. This 
has the effect of introducing an artificial viscosity to the solution proce- 
dure, which assures its stability in supersonic regions where the equations 
are hyperbolic. This artificial viscosity approach is employed principally by 
investigators using the non-conservation form of the full potential equation. 
It has been used in the work of Dulikravich [25 j, Ives [19], Rae [21], and 
Fruhauf [28]. Dodge [18], who also solves the non-conservation form, 
constructs hyperbolic and elliptic operators on a near-characteristic grid 
which is updated during the course of calculation. The artificial compressi- 
bility scheme of Hafez, which achieves the same objectives as the artificial 
viscosity approach, has been employed by all investigators solving the conser- 
vation form of the full potential equation, that is, Caspar [22], Farrell 
[27], Ecer [31], and Hirsch [33]. 

Solution of the Algebraic Equations 

Another distinguishing feature between the various methods for solving the 
full potential equation is the technique which is used to solve the set of 
difference equations resulting from either finite differences, finite volumes 
or finite elements. The first of such techniques, and the most common, is 
some form of relaxation procedure, generally successive line over relaxation 
(SLOR). This is a technique employed in the work of Farrell [27], Dulikravich 
[25, 26], Ives [19, 20], Fruhauf [28], Rae [21], and Hirsch [33], Other 
approaches generally employ a form of direct solver to do an inversion of the 
linearized equations. Gaussian elimination is commonly used here as the non- 
iterative part of an overall iterative solution scheme. Ecer [31, 32], 
Laskaris [29, 30], and Caspar [22, 23] all employ this approach. The final 
method used to solve the difference equations is approximate factorization 
techniques, such as those developed by Holst [34, 35]. ADI methods are a 
particular example of these. Such approaches can be up to an order of magni- 
tude faster than traditional relaxation inversion techniques. However, they 
have not generally been applied as yet to turbomachinery cascade problems, 
except by Hirsch in reference [33]. 

Another approach useful for solving the difference equations even more rapidly 
than the approximate factorization methods is the use of multi-grid techni- 
ques. Such techniques were demonstrated for external flow applications, by 
Jameson [36], who extended the methods initially developed by Brandt [37]. To 
date Hirsch [33] is the only one to apply the multi-grid approach to the solu- 
tion of the full potential equation for turbomachinery. 

Results for Full Potential Equations 

Illustrative results will be presented for some of the latest solutions to the 
full potential equation. These results are for 2D and quasi-3D blade- to-blade 
stream surfaces, as well as complete 3D flow passage analyses. Most of the 
references discussed do an excellent job of predicting transonic flow situa- 
tions on compressor and turbine blade rows with Mach numbers below the 1.3 
level . 

To illustrate, results are shown first in Figure 11 for the two-dimensional 
cascade analysis of Dulikravich [26]. Here flow is computed over a cascade of 
symmetric NACA 0012 airfoils at zero angle of attack and zero stagger. 
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gap-to-chord ratio of 3.6, and inlet Mach number of 0.8. These results agree 
very well with the calculations of Caughey, obtained in private communcation, 
for the identical case. Caughey's calculations were performed using the 
methods of Reference [24]. Peak Mach number in this calculation reached a 
value of 1.31, and the sharpness of the shocks shows the excellent results 
which can be obtained with the full potential approach. Very similar results 
for this case have been obtained by Farrell [27], Ecer [31], and others. 

Results of Farrell [27] are shown in Figure 12 for a supercritical compressor 
stator tip section designed for NASA Lewis by Sanz using a hodograph technique 
based on Bauer, et. al. [38]. The trailing edge of this blade ends in a cusp. 
The inlet Mach number is 0.71, and inlet flow angle 31.16°. The results 
show that the rapid compression is captured very accurately without any over- 
reaction or steepening by the potential method to form a shock. 

Figure 13, again from Farrell, indicates the importance of quasi -3-dimensional 
streamtube effects in the transonic flow regime. This figure shows a series 
of calculations performed on a thick compressor stator hub section, again 
developeo for NASA Lewis by Sanz. Farrell's approach was first used to calcu- 
late strict 20 flow over the blade. This agrees very well with the Sanz hodo- 
graph solution, except at the trailing edge, where in Farrell's analysis the 
idealized blade of infinite length was replaced by a blade with constant 
trailing edge radius. The remaining two curves on the figure show the effects 
of radius change and streamtube convergence. In the first of the remaining 
two curves, a streamtube convergence was imposed of approximately 14% (axial 
velocity density ratio = 1.15), with no radius change on the stream sheet. 

The streamtube convergence strongly increases the Mach number on the suction 
surface of the blade, so that the presence of a reasonably strong shock is 
evident. In the final curve, the opening up of the passage due to a radius 
change of 5%, entirely relieves the effect which was evident previously. The 
increasing radius has a strong decelerating effect on the flow since the 
blade-to-blade passage now diverges in the downstream direction. The differ- 
ences evident in these three curves indicate very strongly the requirement to 
consider quasi-3-dimensional effects in transonic design situations. 

Finally, two results from Hirs'“h [33], shown in Figures 14 and 15, indicate 
typical results obtaineo with a finite element blade-to-blade code. The 
first. Figure 14, shows results calculated for a VKI-LS59 gas turbine cascade, 
calculated on the grid shown in Figure 5. The inlet Mach number is 0.281, 
inlet flow angle 30.0°, outlet Mao'! number 0.975, and outlet flow angle 
-65.89°. The stream channel convergence in the through flow direction is 
unity. The comparison with experimental data in this accelerating flow situa- 
tion is extremely good, even in the transonic regime at the rear end of the 
blade. The compressor results shown in Figure 15 are for a double circular 
arc 9.5° camber compressor cascade. Inlet Mach number in this case is 1.05, 
inlet flow angle 58.0°, outlet Mach number 0.761, and outlet angle 49.5°. 

The stream channel convergence from inlet to outlet has a value of 0.86. The 
deviations in this case on the suction side of the blade between the computa- 
tion and experiment are explained by the presence of large boundary layers 
which exist with compressor flows. 

Finally three-dimensional results computed with the 30 code of Dulikravich 
[26, 39] are presented in Figure 16. These results are for an idealized rotor 
with the flow conaitions illustrated on the figure. The results show regions 
of supersonic flow on all blade surfaces from hub to shroud, with fairly 
strong shocks existing in the tip region. Peak Mach numbers at the tip reach 
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values of approximately 1.5. Two grids were used to obtain this solution, 
progressing from a rather coarse grid with 24 points around the blaoe by 6 
normal and 6 radial to a final grid with double those mesh dimensions in all 
directions. 


STREAM FUNCTION EQUATION ANALYSES 

In the previous section on the full potential equation, it was indicated that 
major work in the application of that equation to internal turbomachinery 
flows didn't occur until after the papers by Murman and Cole [14], Jameson 
[15], and Hafez [16]. These papers appeared in the early and mid 70' s, so 
that all of the major applications to turbomachinery using the full potential 
equation have occurred from the mid 70' s to the present time. A major point 
driving the development of these analyses was the need to extend the 
performance of modern day turbomachinery blading so that significant ranges of 
the flow field were in the transonic regime. 

Stream function equation applications to turbomachinery, on the other hand, 
occured about ten years prior to those for the full potential equation. The 
earliest applications appeared in the late 60' s, and these have continued to 
be developed and extended throughout the 1970's. It may seem strange that 
developments of the stream function equation appeared sooner than those for 
the full potential equation, since stream function analysis permits a more 
general modeling of the full rotational flow field, incorporating all the 
generality contained in the full Euler equations. (This capability is, 
however, limited strictly to subsonic flows, at least up until the present 
time). A major development paving the way for analyses of turbomachinery 
flows using the stream function equation was the series of early papers by Wu, 
particularly [40], which derived the stream function equations for solutions 
on hub-to-shroud and blade-to-bl ade stream surfaces of turbomachinery. This 
paper appeared in 1952, but its application was delayed until the late 60's 
awaiting the development of larger and faster computers on which the equations 
could be accurately solved. 

In Imu's so-called general theory, the stream function equations are solved on 
two intersecting families of stream surfaces - the hub-to-shroud and the 
bl ade-to-blade. The complete three-dimensional flow field is calculated 
through an iterative process which involves transfer of information back and 
forth between these two sets of surfaces. Wu's analysis assumes that the flow 
relative to the blade rows is steady; however it does not demand that the flow 
at the exit of the blade row be axisymmetric. Therefore, exit flows can vary 
circumferentially, and following blade rows are thus subject to time-varying 
inlet flow. Because of this, Wu's general method of analysis is only 
applicable to flow through an isolated blade row. To apply it to a 
multi-stage machine, the time dependence must be removed by circumferentially 
averaging between each row of blades. 

Several principal authors have developed computer codes for quasi-3D 
calculations in turbomachinery, following the work of Wu [40]. Adopting the 
nomenclature of Wu, these flows are either solved along meridional types of 
surfaces, which Wu calls S2, or blade- to-blade surfaces, called SI. Most 
authors have used the approach of having only a single S2 meridional stream 
surface along which the through-flow is calculated, interacted with a series 
of axisyinnetric blade-to-blade surfaces from the hub to the shroud. 
Complications arise in the interaction between these two types of surfaces 
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from the fact that the shape of the S2 meridional stream surface as well as 
its thickness can only really be determined from the SI blade-to-blade 
solutions. These solutions in turn require the knowledge of the through-flow 
on the meridional surface in order to provide boundary conditions as well aS 
the thickness of the blade-to-blade stream sheets. Obviously the procedure to 
approach convergence is iterative, and can be handled in a wide variety of 
different ways. 

Derivation of the stream function equations on the two intersecting surfaces 
is complex, and the reader is referred to details in Wu and in the other 
references to be given shortly. Each author approaches derivation in a 
slightly different way, and the resultant non-conservation equations on the 
meridional surfaces ana blade-to-blade surfaces differ from each other in a 
variety of ways. 

Historical Development and Applications 

The two earliest developers of stream function analyses for turbomachinery 
workeo independently of each other, and their methods appeared at 
approximately the same time. These were Katsanis in the United States and 
Marsh in England. Oddly enough they worked on different surfaces, Katsanis on 
blade-to-blade methods, and Marsh on hub-shroud stream surfaces. Katsanis 
publisheo several NASA TN's for incompressible and high subsonic blade-to-blade 
applications, and then in 1969 published [41] his TN for what he called 
transonic blade-to-blade flows. This analysis applied to any fixed or 
rotating axial, radial or mixed flow turbomachinery blade row. Quasi-3D 
effects were incorporated through a stream channel thickness. The transonic 
flow referred to in Katsanis's title was not obtained by the stream function 
method, but by application of a velocity gradient approach, using information 
about the shape of the streamlines obtained from the stream function method at 
reduced weight flow. Wood [42] has since devised methods for extending 
Katsanis's approaches into the low transonic regime, obtaining more accurate 
results without having to employ the velocity gradient approach. In 1969 
Katsanis and McNally [43] also published a method to analyze blade-to-blade 
flows through slotted or tandem blade rows, as well as a method [44] to 
greatly magnify the solution obtained with the methods of [41, 43] around the 
blade edges or the slot region of a tandem or slotted blade. All of these 
methods of Katsanis use a regular rectangular mesh on the blade-to-blade flow 
surface. Such a mesh has non-uniform-length mesh spacing adjacent to the 
blade surfaces, requiring special treatment of the boundary conditions at 
these locations. 

Marsh [45] was developing at the same time his 20 stream function analysis for 
hub-shroud surfaces. This was the first such analysis to appear in the 
literature; and its major contribution, other than the solution method itself, 
was the development of techniques for applying the method on an irregular 
grid. This grid was composed of parallel lines in the radial direction in 
conjunction with through-flow lines which follow the shape of the hub and 
shroud boundaries. The net result was irregular, or non-rectangular, mesh 
cells in the solution domain. Marsh's technique applied to axial, radial and 
mixed flow turbomachines, and the finite difference equations . for the stream 
function were solved by a matrix method; hence the "matrix through-flow" label 
given to Marsh's techniques and other techniques such as his which followed. 
Marsh deviated from the "general theory" of Wu, as most other authors were 
also to do subsequent to him. Marsh developed what he called the 
"through-flow theory" in which the time dependence was removed by treating the 
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flow on the midchannel as an axially symmetric flow between each pair of 
blades in the blade row. 

Marsh's code was applied at the National Gas Turbine Establishment, and Smith 

[46] presented a paper in 1968 contrasting results from that technique to 
those of the established streamline curvature methods. The major conclusion 
of the paper was that the matrix through-flow method enabled significant 
advances to be made in the calculation of quasi-3D flows. It pointed to the 
need for incorporating a good endwall boundary layer solution and accurate 
loss mechanisms in this sort or program before adequate comparisons could be 
made with experimental data. 

Smith and his co-workers at N.G.T.E. continued the development of matrix 
methods and quickly extended them to blade-to-blade surfaces. In 1970 Smith 

[47] and Frost presented a paper for computing flow fields on blade-to-blade 
stream surfaces using both a matrix stream function analysis and a streamline 
curvature technique. This method was applicable to any type of axial or mixed 
flow compressor or turbine blading with either stationery or rotating blade 
rows. In 1970 Smith [48] also published a paper describing both the 
meridional and the blade-to-blade analyses which were in use at that time at 
N.G.T.E. This was the first paper describing the meridional and blade-to-blade 
methods being used together in a unified way. However, the paper did not 
describe in detail the iteration process used between the two approaches. 

Another set of codes for both hub-to-shroud and blade-to-blade analyses using 
the stream function equation was developed in the early 70's at Carleton 
University in Canada by Davis and Millar. Intially Davis [49] published a 
thorough analysis of Marsh's approach to generating a curvilinear grid on the 
hub-shroud stream surface. His paper redeveloped the Marsh technique and 
again applied it to the hub-shroud stream surface of an axial turbomachine. 

In 1972 and 1973 Davis and Millar [50, 51] published extensive reports on the 
development of both blade-to-blade and hub-to-shroud codes at Carleton 
University. The hub-shroud code, extended the methods originally developed by 
Marsh. The blade-to-blade analysis was likewise similar to those published 
earlier by Katsanis and Smith. In 1974 Davis and Millar [52] compared the 
matrix through-flow technique to streamline curvature methods for calculating 
flows on hub-shroud surfaces. They applied the two techniques to a duct flow, 
a transonic fan, and three-stage axial compressor. These comparisons 
indicated that the two approaches gave similar results, and that there was a 
small operational aovantage with the matrix through-flow method. Both methods 
were shown to be subject to certain instabilities. The convergence of the 
stream function from iteration to iteration had to be damped in the matrix 
through-flow method, while the shift in streamline position had to be carefully 
damped in the streamline curvature techniques. 

In 1975 Davis [53] presented the first hub-to-shroud stream function solution 
for flow in a centrifugal compressor. To accomplish this he used a special 
form of through-flow grid which remained quasi-orthogonal throughout the 
solution domain. He employed two versions of the stream function equation, 
one for regions where the axial velocity exceeded the radial, and the second 
in regions were the radial velocity was the larger. Furthermore, Davis 
incorporated a turbulent endwall calculation, using an integral method based 
on the entrainment theory of Head [54]. Although the equations derived by 
Davis were for application to centrifugal turbomachinery, they were only 
applieo in the paper to stationary components: an inlet, a diffuser, and an 

intra-stage return bend. 
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In 1973 and 1974 Katsanis and McNally published NASA reports and an 
accompanying paper [55] describing their hub-to-shroud stream function 
analysis for use with the blade row analysis of reference [41]. This 
hub-shroud analysis could handle either axial or mixed flow compressor or 
turbine blade rows, but was limited to a single blade row. Later, in 1977, 
the same authors [56, 57] extended this code so that it also applied to radial 
or centrifugal blade rows. Mild transonic flows were treated with the same 
velocity gradient method used in the blade-to-blade analysis. 

At about this same time. Marsh [58] published a paper similar to earlier 
papers by Smith [46] and Davis [52], comparing the matrix through-flow 
analysis approach to streamline curvature techniques for the hub-shroud 
problem. Like the previous papers, he discussed various loss models which 
were in use, and the need for a good endwall boundary layer analysis. Marsh 
concluded that the matrix through-flow and the streamline curvature techniques 
could be viewed as two different methods for solving the same governing 
equations on the same mean stream surface. He did not conclude that there was 
a definite superiority of one method over the other. He recommended that work 
be pursued to develop more accurate methods for estimating the losses within 
blade rows, to calculate the development of the endwall boundary layers, and 
to predict secondary flows. 

Iterative Approaches and Finite Element Methods 

In 1976 the first paper was published, by Bosman [59], giving significant 
detail concerning an iterative approach to couple hub-to-shroud and 
blaoe-to-bl ade stream function analyses. Bosman presents equations applicable 
to any type of stationery or rotating turbomachinery geometry including 
centrifugals. He described the following iterative procedure: The initial S2 

stream surface is assumed to coincide with a mean blade shape. From the 
calculated S2 streamlines, SI stream sheets are generated by revolution about 
the axis of rotation. SI solutions are then obtained and from these 
mass-averaged streamlines are defined, thus giving a re-definition of the S2 
stream surface. Another S2 solution is then obtained, and the process 
repeated until convergence is obtained. Thus the shape of the S2 stream 
surface evolves from SI solutions and is not constructed by any sort of 
corrections to the mean camber surface for incidence and deviation angles as 
is done in several other techniques. Bosman applied his techniques to 
calculate flows in a low speed centrifugal compressor, ana in a radial inflow 
turbine. A meridional view showing the hub and shroud profile for the low 
speed compressor as well as the grid used in Boman's calculation is indicated 
in Figure 17. This is an eight bladed compressor with a very short inducer 
section. In Figure 18 the SI mean, streaml i nes calculated by Bosman at the hub, 
mid-section and shroud for approximately zero incidence are shown. It is 
unlikely that a user trying to predict the stream surface to satisfy given 
incidence and slip conditions would have produced streamlines at all 
approximating those calculated by Bosman. Figures 19 and 20 show the 
calculated hub and shroud S2 velocity distributions, and illustrate the effect 
of the stream surface shape upon these. These figures show that there is a 
major influence on the hub and shroud velocities due to the graduated slip 
that occurs in the exit sections of the impeller. 

A second, more elaborate iterative procedure was described by Adler and 
Krimerman [60] in 1978. Their meridional and blade-to-blade computer codes, 
based on finite element methods, are described in references [61, 62]. 

Adler's iterative process occurs in four distinct steps. The first two are 
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normal applications of meridional and blade- to-blade calculations on 
axisymmetric stream surfaces. In the initial hub-shroud solution the mean 
camber surface is adjusted at the loading edge to be aligned to the inlet 
relative flow velocity and at the trailing edge to be parallel to an assumed 
deviation angle. After these initial two solutions a number of different 
hub-to-shroud stream surfaces are calculated from corresponding streamlines in 
the blade-to-blade surfaces. These hub-shroua surfaces are no longer 
axisymmetric. From these multiple hub-shroud solutions, blade-to-blade 
surfaces are finally calculated by connecting corresponding streamlines. 

These blade-to-blaae surfaces are therefore no longer surfaces of revolution. 
Iteration is continued back and forth between the third and forth steps until 
convergence is obtained. Adler applies this technique to a centrifugal 
impeller, and shows that the results clearly deviate from results obtained 
using a single axisymmetric meridional surface. He claims that these 
deviations are significant enough to justify this type of flow field 
calculation in highly-loaded compressors with back-swept blading where the 
flow is very three- dimensional. Since the four corner streamlines in Adler's 
analysis are forced to remain within the corners, one can argue whether the 
results obtained here are worth the extra effort. 

Hirsch [63] in 1976 published the first numerical solution of the meridional 
through-flow stream function equation based on the finite element method. The 
method was applied to axial flow machines, but was also developed for radial 
machinery. The method was shown to be applicable to transonic stages in cases 
where the tangential velocity distribution was given, as long as the meridional 
velocity remained subsonic. Several years later, Hirsch published an iterated 
analysis [64] in which meridional and blade-to-blade finite element stream 
function analyses were combined for application to axial turbomachinery. 

Again the flow in the S2 surface was replaced by the calculation of the exact 
mass-averaged pitch-averaged flow on the meridional ( ) plane. Hirsch 

uses second-order isoparametric quadri laterial elements which allow an 
accurate simulation of blade curvature even in highly curved regions such as 
leading and trailing edges. Even though Hirsch's meridional calculation 
permits transonic relative flows as long as the merdional Mach number is lower 
than 1, the blade-to-blade code used in this combineo analysir permits only 
subsonic velocities throughout the flow field. Results of this analysis were 
compared favorably with the LUV data obtained by DFVLR in Germany [65]. In 
1980 Hirsch [66] again presented a combined iterative analysis, this time for 
centrifugal compressors. Results were computed for the radial compressor 
mapped with LDV by Eckardt [67]. Hirsch concluded that the viscous and 
secondary flow effects were not well reproduced with the through-flow 
calculations particularly in the back ena of the compressor. Hirsch also 
presented a quasi-3D calculation on the so-called Type-B centrifugal 
compressor described in reference [68]. Figure 21 shows the shape of this 
compressor in the meridional plane and indicates the finite element mesh as 
well as the calculated meridional streamlines obtained for a flow coefficient 
of 0.2. Figures 22 and 23 show the static pressure distributions on the hub 
and shroud sections for a flow coefficient of 0.5. The experimental results 
from [68] are compared here to the calculated blade pressure distributions. 

The inviscid calculations predict a stronger local acceleration along the 
shroud suction surface than experimentally determined. Figures 24 and 25 show 
the calculated Mach number distributions along the hub and shroud. The figure 
at the hub indicates a local low velocity region in the outlet portion of the 
blade near the pressure side of the passage. The shroud profiles indicate 
details ot a local acceleration region founo near the front end of the suction 
surface of the blade. The differences between data and calculation indicate 
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that the secondary flows existing in the passage create velocity components 
which give a better behavior to the real flow than that found in the 
calculations. These differences illustrate the limitations of the quasi-3D 
approach were secondary flows are not taken into account. 

Finally, in 1980 Goulas [69] presented a stream function analysis for the 
blade-to-blade flow in a centrifugal compressor which contains splitter 
blades. He derived his analysis so it could either calculate isentropically, 
or simulate turbulent flows with the addition of a simple turbulence model and 
zero velocity conditions at the wall. He adapted the method to handle 
stagnation points as well as the formation of small re-circulation zones at 
the front end of the splitter blades which sometimes occur in the analysis. 

His method was applied to a centrifugal compressor in which he studied various 
axial locations for the leading edge of the splitter. 

Multi-Stage Meridional Capability 

Of the methods just described for analyzing hub-shroud flow with stream 
function analyses, the methods of Marsh [45], Smith [48], Davis [49, 51], and 
Hirsch [63, 64] all permit the analysis of multi-stage machines. The method 
of Davis, however, was only demonstrated [49, 51] for a single full stage 
machine. On the other hand, the methods of Katsanis [55, 56], and Bosman [59] 
only apply to a single blade row. 

Generation of Algebraic Equations 

Most of the methods described above for both meridional and blaoe-to-blade 
types of stream function analyses have used finite difference techniques to 
generate the algebraic equations from the partial differential equation. The 
methods of Marsh, Smith, Davis, Bosman, Katsanis, and Hafez all fall within 
this class. Only the more recent methods of Hirsch and Adler use finite 
element techniques to discretize the equations. As time goes by and more and 
more automatic grid generators are developed for finite elements in 
turbomachinery applications, this approach will no doubt be used more heavily. 

Iteration Between Hub-Shroud and B1 ade-to-Bl ade Solutions 


Most of the authors referenced above, except Marsh who developed the first 
through-flow method, have developed codes for both meridional and 
blade-to-blade analyses. In some of these references details were given 
concerning automated interactions between the SI and S2 surfaces. Marsh, 
Smith, Davis, and Katsanis described no such interactions. Bosman [59], Adler 
[60], and Hirsch [64, 66] all described such an interaction between the two 
types of surfaces. In cases where no interaction occurs, the midchannel 
stream surface is either formed as a midchannel projection of the mean camber 
line of the blade, or by altering that projection to accommodate at the 
leading edge for incidence and at the trailing edge for deviation angle. 

Where interaction does occur, the shape of the meridional stream surface is 
obtained from some sort of integration or mass averaging of the flow on the 
blade-to-blade surfaces from hub to shroud. 

Type of Finite Difference Mesh 

In all of the hub-to-shroud analyses mentioned above, there is a fair degree 
of uniformity in the grids, with most using a form of quasi-orthogonal grid 
similar to those originally chosen by Marsh in the first meridional analysis. 
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None are rectangular, and all have some sort of through-flow streamlines which 
match the hub and shroud geometries coupled with a grid of lines nearly 
orthogonal to these passing from hub to shroud. Generally details are not 
presented concerning the generation schemes used for these grids, except by 
Katsanis who uses the method of reference [70]. In some of the later methods, 
particularly those due to Bosnian, Adler and Hirsch, the grids coincide with 
the blade edges giving a more accurate solution in those regions. 

On the bl ade-to-bl ade surfaces a wider variety of meshes are used. In the 
blade-to-blade methods of Katsanis [41], and bosman [59] rectangular grids are 
used with the inherent problems these bring near the boundaries. Rectangular 
grids provide a great advantage in the interior regions of the flow domain. 
Because of their regularity, the finite difference expressions are very simple 
and uniform throughout such a grid. However, at the boundaries such grids 
pose difficulties because of the irregularly sized mesh legs which intercept 
the boundaries in both directions. In other methods, noteably Smith's [47,48] 
and Davis [50], a channel or H-type grid is used with straight lines running 
in the 6 direction and through-flow streamlines adhering to the shape of the 
blade surfaces on the suction and pressure sides of the channel. Finally, for 
the finite element methods, Hirsch uses second-order isoparametric 
quaori laterial elements in both his meridional and blade-to-blade solutions. 
Adler, on the other hand, uses linear triangular finite elements. 

Solution of the Algebraic Equations 

Two principal techniques were used to solve the algebraic equations in most of 
the references. The first involved classical matrix inversion techniques and 
the second relaxation procedures. Of those that mentioned their methods for 
solving the equations, only Katsanis used a relaxation method. Most of the 
authors using direct inversion factored their matrices to lower and upper 
triangular banded matrices. Once this is done, and the banded matrices are 
stored, they can be applied efficiently from iteration to iteration without 
any large expenditure of computer time. 

Viscous Loss Corrections 


A nuniber of different approaches are used with both the meridional and 
blade-to-blade stream function analyses to simulate the effects of viscous 
loss in these inviscid calculations. The first of these involves correlation 
for loss as a function of various blade geometry, setting angle and loading 
parameters. The second method involves a pre-specification of a distribution 
of total pressure from inlet to outlet through the blade row. This total 
pressure is incorporated into the solution procedure and effects the velocity 
accordingly. The final method employs the calculation of either endwall or 
blade surface boundary layers. These boundary layers then alter the shape of 
either the meridional or blade-to-blade channel thus effecting the internal 
flow solution. 

Transonic and Three-Dimensional Flows 

The stream function method has always traditionally been applied to 
two-dimensional and subsonic flows. Recently Hafez [13] has approached the 
solution of the stream function equation in conservation form using techniques 
traditionally applied to the full potential equation in that form. He has 
found that traditional methods such as artificial compressibi 1 ity can be 
extended to the stream function equation. He has investigated methods to 
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overcome the fact that the density is not uniquely determined in terms of the 
mass flux, having instead both subsonic and supersonic solutions on either 
side of the sonic point. He has applied the resultant method to flows over a 
NACA 0012 airfoil at a variety of inlet Mach numbers in the transonic range. 

He has also performed calculations about a 10% parabolic airfoil, as well as 
transonic flows about a cylinder. These results have been compared to both 
potential and Euler solutions. Hafez has likewise investigated the 
application of stream function methods to three-dimensional flows through the 
use of two different stream functions for this problem. This work is ongoing, 
and as yet has not been applied to turbomachinery situations. 


FULL EULER EQUATION ANALYSES 
Explicit Time Marching Methods 

Explicit solutions to the Euler equations have been under development for a 
number of years, with applications for internal flow situations dating back to 
the 1950's. An explicit method is one in which all spatial derivatives are 
evaluated using known conditions at an old time level. Futhermore, 
information at the new time level depends on information obtained from only a 
small number of points. The resultant methods are simple and easy to code. 

All such methods, however, are limited by the so-called Courant, Friedrichs, 
and Lewy (CFL) stability limit, which states that the domain of dependence of 
the numerical finite oifference scheme must contain the complete domain of 
dependence of the original hyperbolic differential equations. 

A major milestone in the development of explicit methods was the paper by 
MacCormack [71]. This method has second order accuracy in both time and 
space. It is a two-step predictor-corrector method, which alternates between 
forward and backward differencing on the two steps. The ease with which this 
method can be applied has led to many applications in the turbomachinery 
field. One of the first applications of the MacCormack scheme to 
turbomachinery was by Gopalakri shnan and Bozzola [72]. Gopalakrishnan applied 
the basic MacCormack algorithm without modification to a transonic compressor 
cascade with supersonic inlet flow shocking down in a blade row to subsonic 
outlet flow. 

Another application of the same scheme to turbomachinery was that of Kurzrock 
and Novick [73]. This solution was obtained for a rotating blade-to-blade 
stream surface, with radius change and stream channel convergence included. 

The authors retained the viscous terms from the Navier-Stokes equations, with 
an artifically enlarged viscosity coefficient in order to capture shocks. 

They applied the method to transonic flow in a 20 compressor cascade and also 
to the quasi-3D stream surface of a compressor rotor. Comparisons between 
computed and experimental exit conditions were presented. 

Thonipkins [74] has applied the MacCormack algorithm to flow through a three- 
dimensional transonic compressor rotor. This method can be applied to any 
general compressor blade shape, including those with part-span shrouds. 
Computer results applying this method to a transonic compressor will be 
presentee later. 

Another set of explicit methods is obtained by writing the conservation laws 
in integral form and applying them to local control volumes surrounding each 
grid point. The fluxes of mass, momentum and energy crossing the control 
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surfaces are evaluated using the values from surrounding points. A first 
order integration in time is then used to advance the dependent variables 
forward to a steady state. 

The first major application of this technique for turbomachinery was by 
McDonald [75], who applied the method to 2D transonic flow in axial turbine 
cascades. His approach includes the use of hexagonal elements surrounding 
each grid point, as well as the replacement of the energy equation by a 
statement of constant total enthalpy throughout the field. Computed and 
experimental results were compared for a number of high-turning turbine 
cascades with supersonic exits. 

Denton [76, 77, 78] has developed a somewhat simpler method for both 2D and 3D 
turbomachinery flows. Denton employs quadrilaterial elements in two 
dimensions and six-sided elements in 3U, which lead to simpler expressions for 
his surface flux integrals. In order to ensure stability he uses upwind 
differencing in the streamwise direction for the fluxes of mass and momentum, 
while using downwind differencing for pressure. Centrfal differencing is used 
for all quantities in the pitchwise direction. This scheme has the property 
that stability depends only on the axial Mach number, not on the absolute Mach 
number which is more usual. This method has been successfully applied to a 
wide variety of both axial and mixed flow compressor and turbine geometries. 

Another first-order explicit method has been developed recently by bosman and 
Highton [79] for three-dimensional flow situations. The method employs two 
separate overlapping grids on which density and internal energy are evaluated 
at one set of nodes, and velocities are evaluated at the second set. This 
overlapping scheme facilitates the evaluation of fluxes for the corresponding 
control surfaces, bosman updates his primative variables in a sequential 
fashion. First, he evaluates the velocities, then density and pressure. A 
new set of velocities is then obtained using the new pressures. Next, 
internal energy and temperature are updated, which leads to a final update of 
the pressure. At each stage in this process, the most recent values of all 
primative variables are used to update the fluxes. The method has been 
applied to both radial inflow turbines [79] ana to centrifugal compressor 
impellers [80]. Results from the radial turbine example will be presented 
later. 

Recent efforts have been devoted to improving both the accuracy and the speed 
of explicit methods. A significant example of the former is Moretti's X - 
scheme [81] which exploits concepts from the method of characteristics for 
hyperbolic systems. In Moretti’s scheme he rewrites the Euler equation so 
that the right-hano sides involve derivatives of one-dimensional Riemann 
invariants only. These derivatives are then replaced by one-sided differences 
in directions corresponding to the projection of the associated 
bicharacteristics onto the previous time plane. The equations are then 
updated in time by a two-step predictor-corrector scheme similar to that of 
MacCormack ' s . The result seems to be improved accuracy as evidenced by very 
sharp shocks captured with mooest numbers of grid points. The A -scheme by 
itself only produces isentropic shocks unless it is corrected with a shock 
fitting procedure, as is done in De Neef and Moretti [82]. The X -scheme 
has been modified and applied to simple compressor and turbine cascades by 
Pandolfi and Zannetti [83]. A similar method to the X -scheme has been 
developed by Chakravarthy, et al. [84]. 
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Ni [85] has developed a new explicit method which is both accurate and quite 
fast. He begins with a scheme that appears to be equivalent to the second- 
order Lax-Wendroff procedure, (see Richtmyer and Morton [86]). However, by 
using a spatially varying time-step, which is taken to be everywhere near the 
CFL limit, he obtains a method which operates with the largest possible time 
step in all regions. More importantly, this effectively biases the 
differencing in such a way that the finite difference scheme has a domain of 
dependence which approximates that of the underlying hyperbolic system. Ni's 
method has been coupled with a multi-grid procedure which greatly speeds 
convergence to a final steady state. The method was applied in [85] to 
transonic flow in a turbine cascade as well as to an axisymmetric nacelle with 
centerbody . 

Another approach to speeding convergence of the explicit methods to a steady 
state is to add purely artifical unsteady terms to the steady equations. When 
properly constructed, these artificial terms can introduce a strong internal 
damping into the resulting unsteady system. Such approaches are called 
pseudo-unsteady. One such method has been employed by Essers [87] to compute 
2D steady irrotational transonic flows. In this approach the artificial 
unsteady terms are added to the continuity equation and to the irrotational ity 
condition. This results in two equations in the unknown velocities, with 
density obtained from the isentropic relationship. These equations are solved 
by an adaption of MacCormack's predictor-corrector scheme. Results are 
presented by Essers for two-dimensional flow through a turbine rotor blade 
section with supersonic exit. Effects of various treatments of the blunt 
trailing eoge are also presented. 

Another pseudo-unsteady method has been developed by Viviand and Veuillot [88, 
89]. In this technique the energy equation is replaced by the statement of 
constant total enthalpy, and the pressure is then expressed as a simple 
function of density and velocity. The resultant system includes the unsteady 
continuity and momentum equations. These are solved by a generalization of 
MacCormack's predictor-corrector scheme. Since the system has no unusual 
damping mechanism beyond the normal damping due to truncation errors and 
simple artificial viscosity, it relies on careful treatment of waves at the 
boundaries in order to reach a steady state. A three-dimensional version of 
this method has been developed by Brochet [90, 91] and applied to flow in a 
supersonic compressor cascade with converging endwalls and to transonic flow 
in a fan rotor. Results from the cascade will be presented later. 

Implicit Time Marching Methods 

Implicit time-marching methods for both the Euler and the Navier-Stokes 
equations date from the mid-1970's. In implicit methods the equations are 
backward differenced in time, and the non-linear terms at the new time are 
Taylor-expanded about their values at the previous time level. This produces 
a system which is linear in the unknowns at the new time level. The spatial 
derivatives are then approximated by finite differences, resulting in a large 
system of coupled linear algebraic equations for the unknowns at the new time 
level. In each of the methods to be discussed here, these equations are 
solved by block alternating-direction-implicit (block ADI) techniques. In 
this approach the matrices are first factored into a sequence of matrices for 
one-dimensional problems, each of which can be inverted by a tri-diagonal 
routine. The matrix elements in these one-dimensional problems are in turn 
simple block matrices whose size is equal to the number of unknowns at each 
grid point. The solution at each time-step proceeds non-iterati vely by first 
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moving along the grid lines in one direction, inverting a tri-diagonal matrix 
for each line. It then proceeds similarly/ in each of the other physical grid 
directions of the problem. 

The first of these methods was introduced by Briley and McDonald [92, 93], 
primarily for the compressible Navier-Stokes equations. Beam and Warming [94] 
independently developed a very similar method for the Euler equations. Briley 
and McDonald [95] have since shown that when the Beam and Warming algorithm is 
written in the "delta" form to solve for the corrections to the unknowns at 
the new time levels, the two methods have identical linearized block implicit 
matrices. 

Steger [96, 97] has developed a curvilinear coordinate version of the 
Beam-Warming algorithm for viscous as well as inviscid flows, and has applied 
it to both isolated airfoils and cascades in two dimensions. Shamroth, et al . 
[98] has applied the Bri ley-McDonald procedure to laminar and turbulent flow 
through a cascade. Finally, Fruhauf [28] has applied the Beam-Warming 
algorithm to solve the Euler equations for both subsonic and supercritical 
flow through cascades. 

In all of the time-marching explicit and implicit methods some form of 
numerical damping is present in the solution procedure which smoothes the 
oscillations that occur in the vicinity of strong shock waves. Most of the 
methods have added the damping explicitly in the form of a higher-order 
differencing term. However, in all of these cases the natural truncation 
errors that occur as the result of any finite difference procedure will add 
their own damping. Since the form of these damping terms varies considerably 
from method to method we will not discuss the details for any particular 
appl ication. 

New Methods Under Development 


A number of methods are under development in order to achieve more accurate 
and faster solutions for the Euler equations. Delaney [99] is developing a 
hopscotch method for solving the Euler equations for application to cascades. 
The method appears to be significantly faster than the original MacCormack 
algorithm. 

Denton [lOO] has extended his earlier Euler method by employing a simpler more 
accurate differencing scheme. He has also increased the convergence speed 
through the use of a simple multi-grid procedure. He has applied the method 
to 2D transonic flow in both compressor and turbine blade rows. 

Johnson [101], in order obtain the benefits of existing rapid solution 
procedures for second-order partial differential equations, has developed a 
new technique in which the first-order steady Euler equations are imbedded in 
a second-order system. Published results to date have been for the transonic 
small disturbance equations and the full Euler equations in subcritical flow. 

Both Ecer and Akay [102] and Lacor and Hirsch [103] have developed methods 
which are the equivalent of solving the full steady Euler equations. In these 
two approaches the velocity is split into potential and rotational parts. The 
resultant system of equations includes a second order equation for a potential 
function which is obtained from the continuity equation, and a pair of first 
order convective equations describing the evolution of two scalars which 
together determine the rotational part of the velocity field. These are 
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solved as a coupled system. Both of these methods have used finite element 
formulations. Ecer has applied his technique to 20 transonic flow in a 
channel with a bump and to flow around a 20 cylinder. Hirsch applies his 
method to calculate 30 flow in a rectangular elbow with 90° of turning. 

Finally Chang and Adamczyk [104] have developed a new semi-direct algorithm 
for computing three-dimensional inviscid shear flows. This algorithm is 
composed of two iteration loops. In the inner loop, the velocity and density 
fields are evaluated for specified vorticity, total enthalpy and entropy 
fields. This evaluation is reduced to the solution of a pair of Poisson 
equations in the computational domain which are solved by three-dimensional, 
fast, direct Poisson solvers. In the outer loop, the vorticity, total 
enthalpy and entropy are obtained by solving convective equations for a pair 
of scalars in a manner similar to that of the previous two references. In the 
present work, finite difference procedures are used throughout. To date the 
method has been applied to study the development of inviscid shear flows in 
turning channels. 

Applications of Euler Methods 

Several results will be presented to indicate the kinds of calculations which 
are being performed with the various Euler equation methods. The method of 
Thompkins [74] has been applied at NASA by Chima and Strazisar [105] to 
calculate the three-dimensional flow field within a transonic axial compressor 
rotor at design speed, and to compare those results to laser anemometer 
measurements at maximum-flow and near-stall operating points. Figure 26 
indicates Mach number contours for the measured laser anemometer results at a 
section 15% from the tip of the blade. These results can be compared to the 
calculated contours in Figure 27 at the same location. These figures indicate 
a pronounced bowwave and passage shock system, and show excellent agreement 
between the measured and calculated results. 

The method of Bosman [79] has been applied to a radial inflow turbine which 
turns through 90° of deflection in the meridional plane and has a 70° 
deflection in the vvi- G plane at the blade tip. Such a geometry will 
naturally produce large three-dimensional secondary flows, and one of the 
strong points of the Bosman tecimique is that it picks up these natural 
inviscid vortex motions in such a three-dimensional geometry. Figure 28 
indicates for the design point condition the calculated streakline pattern on 
the hub section of the blade. Some of these streaklines have been joined to 
form a streamline which is seen to migrate from mid-passage over to the 
suction surface of the blade. A migration in the opposite direction was 
calculated by Bosman along the tip. Similar calculated patterns are indicated 
for the suction ana pressure sides of this turbine rotor in Figure 29. Figure 
30 compares the mid-passage-flow shroud static pressure distribution with 
experimental results and with two-dimensional calculations using a blade-like 
hub-shroud stream surface. The two-dimensional results have strong deviation 
from the experiments, especially in the trailing edge region, while the 
three-dimensional calculation agrees with the experiments quite well. 

Brochet [90] applies his method to the calculation of flow through a 
supersonic compressor cascade with subsonic axial velocity. This cascade has 
converging sidewalls in the throughflow direction, with a contraction ratio of 
0.7. The calculation was obtained for an upstream Mach number of 1.5 and a 
compression ratio of 2.0. Calculated pressure contours, related to upstream 
stagnation pressure, are presented in Figures 31 and 32 for the plane of 
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symmetry section and for the converging wall blade section. For the back 
pressure chosen, a shock is calculated which spans the entire passage from one 
converging wall to the other. At mid-pas. age the calculated flow field and 
shock structure agree quite well with two-dimensional calculations obtained by 
taking into account the stream channel convergence. This is no longer true, 
however, near the converging wall where only the truely three-dimensional 
model appears capable of representing phenomena correctly. 

Finally, the method of Ni [85] has been applied to flow past a VKI gas turbine 

rotor blade, and compared to data by Sieverding [106]. Figure 33 presents 

calculated and experimental surface Mach numbers for this turbine rotor. The 
agreement is excellent on both blade surfaces except for small deviations very 

close to the shock impingement point on the suction side of the blade. Figure 

34 indicates the mesh used in the calculation and shows the calculated Mach 
number contours. It also indicates the experimentally determined shock 
locations, showing that the shock is generated at the trailing edge. 

Calculated and experimental shock locations are in good agreement. 


VISCOUS METHODS 


Full Viscous Equations 


Most flows of engineering interest are adequately described by the 
compressible Navier-Stokes equations. These equations, in two-dimensional 
conservation form and Cartesian (x, y) coordinates, are written 
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where U, F, and 6 are the same as for the Euler equations, and 
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Here JX is the viscosity coefficient and 
conductivity. 



K is 



the coefficient of thermal 


The above equations are valid for turbulent flows, but such computations are 
impractical today due to the large range of length scales in the turbulent 
spectra. The above equations are replaced by time-averaged equations and the 
Reynolds stresses, e.g. are modelled through the addition of 

auxiliary algebraic or differential relations. If a Boussinesq or "eddy 
viscosity" model is adopted for the Reynolds stresses, the above equations 
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will hold for the mean fluid variables f , U. , V" , p , with AX and 
replaced by their effective turbulent values. ' 

For practical engineering flows with curved boundaries, the above Cartesian 
equations are transformed to curvilinear body-fitted coordinates. This adds 
considerable complexity to the equations, especially so for the viscous terms. 

For steady flows the time derivatives are often retained as an aid in the 
solution process. The resulting time-averaged Navier-Stokes equations in 
curvilinear coordinates are quite difficult to solve for a number of reasons. 
First, there are many disparate length scales which must be resolved. These 
may be associated with boundary layers, wakes, vortices and shock waves. For 
complex flows this requires a very large number of grid points even for a 
minimum description of important phenomena. Second, a large number of 
quantities are needed at every grid point in the field. For example, in 
three-dimensional flow, one might need five primative variables, two 
turbulence properties, and nine or more metric derivatives. The above 
requirements lead to the need for a very large computer memory and associated 
long running times to achieve a solution. If sufficient computer core is not 
available, the metric derivatives may have to be re-computed at every time 
step or iteration, which adds to the overall run time. Third, the numerical 
problems associated with the solution of the first order inviscid Euler 
equations are still present in viscous problems especially for high Reynolds 
numbers. One such problem often occurs when the local cell Reynolds number, 

» exceeds 2 and the solution becomes either unstable or highly 
inaccurat^. The common "fix" of locally switching from central to first order 
upwind differencing of the convective term can introduce excessive numerical 
diffusion. Finally, we note in some cases with steady boundary conditions a 
steady solution may not even exist. 

Many techniques have been developed over the years for viscous problems which 
solve simplier sets of equations. We consider only those for steady flows, 
and group them into two categories. 

Partially Parabolic Approximation 

The first of these is the so-calleo partially parabolic approximation. This 
assumes the existence of a predominant flow direction, which is known a priori. 
Therefore, flow separation is excluded. The viscous terms are simplified by 
neglecting diffusion in the main flow direction. This is much like the 
boundary layer approximation. In two dimensions, if x is the flow 
direction, the vector R is dropped from Eq. (20) and d^J^Xis dropped 
from Both terms in & are retained as they are comparable in 

magnitude for compressible flow. For three dimensions see Caretto, Curr, and 
Spalding [107]. One variable, usually the pressure, is treated as elliptic 
and stored at every grid point. The remaining variables, eg. \X, O' , in 20, 
are treated as parabolic and stored on only two or three cross-sections at a 
time. Starting from an approximate pressure field the momentum and energy 
equations are marched in the flow direction. The variables are corrected 
locally to satisfy continuity. After each marching sweep, the pressure field 
is updated by solving an elliptic (Poisson) equation on the entire grid. This 
sequence of marching followed by the pressure update procedure is repeated 
until convergence. This type of analysis is applicable to many internal flows 
in the absence of streamwise separation. This includes those in turning ducts 
and turbomachinery blade rows provioed that the details of the flow in the 
leading and trailing edge regions can be glossed over. The efficiency of such 
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coniputations should compare quite favorably with those of time marching 
procedures if convergence can be obtained in a moderate number of marching 
sweeps, perhaps under 50 for a three-dimensional problem. This would be 
especially true if some rapid procedure such as multi-grid is used to solve 
the elliptic equation for the pressure. 

Fully Parabolic Approximation 

Next we consiaer the fully parabolic approximation. This uses the same 
assumptions as the above with regard to the predominant flow direction and the 
neglect of viscous diffusion in this direction. In addition, upstream 
transmission of pressure disturbances generated during the calculation is 
assumed to be negligible. An initial pressure field, which is stored on the 
full grid, is assumed to contain all of the effects of boundary curvature. 

The remaining variables, including the pressure correction, are stored on only 
two or three cross sections at a time. The momentum and energy equations, or 
an equivalent set, are marched a single time in the flow direction. The 
variables are corrected locally at each cross section in order to satisfy 
continuity. This type of procedure should be applicable to flows in duct-like 
geometries with moderate turning in the absence of streamwise separation. 

Since only a single marching sweep in employed, these computations should be 
orders of magnitude faster than for any time marching procedure. 

V ( 

Both of these parabolic a^iproximations are capable in principle of predicting 
strong secondary flows prbvided that local continuity is well satisfieo on 
each cross section. Both of these can even treat modest amounts of streamwise 
separation if the Flugge-Lotz and Rehyner approximation is adopted, i.e. 
streamwise velocity is artificially prevented from going nega’tive in the 
convective term only. ^ 

We now proceed to discuss methods in each of the above categories starting 
with the fully parabolic .approximation. Within each category, we consider the 
main elements of a solution procedure and in so doing present what in our view 
are the significant distinguishing features of each method. 


FULLY PARABOLIC METHOUS ( 

Main Parabolizing Assumption 

All of the methods in this category require an additional assumption, beyond 
the neglect of streamwise viscous diffusion, in order to obtain a fully 
parabolic system of equations. In the usual procedure a bulk pressure 
correction, , assumed uniform over each cross section, is introduced into 
the streamwise momentum equations. This correction is determined so as to 
ensure the correct mass flux through each cross section. The cross flow 
equations, on the other hand, retain a separate pressure correction which 
is permitted to vary over the cross section. This procedure is employed by 
Patankar and Spalding [108], Briley [109], Ghia et. al . [110], Roberts and 
Forester [111], and Briley and McDonald [112]. 

A different procedure is employed by Anderson [113] in 2D and Anderson and 
Hankins [114] in 3D. The equations are parabolized by writing them in an 
approximate intrinsic coordinate system which could be obtained for example 
from an incompressible potential flow solution. The assumption of small 
velocities normal to the streamwise grid lines eliminates convective 
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derivatives and viscous terms in the transverse momentum equations. The 
resulting system is fully parabolic with characteristic surfaces coincident 
with the cross-planes. Hence, no bulk pressure correction is required. 

Satisfaction of Local Continuity 

We next consider the technique for satisfying local continuity over a 
cross-section. Patankar and Spalding [108] introduce approximate relations 
between the velocities and pressure corrections which are obtained from the 
transverse momentum equations. Substitution of these into continuity yields a 
2D elliptic equation which is solved over the cross-section. Another 
approach, adopted by Briley [109] and Ghia et. al. [110], assumes a 2D 
potential for the transverse velocity corrections. Substitution of this 
potential into continuity yields a 2D elliptic equation for in the cross- 
plane. The divergence of the transverse momentum equations provides a second 
elliptic equation for the pressure correction. Roberts and Forester [111] 
work directly with the divergence of the transverse momentum equations which 
gives a 2D elliptic equation with a source term related to the non-satisfaction 
of local continuity. This equation is solved iteratively with the momentum 
equations until continuity is also satisfied. This technique is related to 
that of Harlow and Welch [115] for 2D time-dependent flows. Briley and 
McDonald [112] split the transverse velocity into irrotational and rotational 
parts described by a 2D potential and stream function .respectively. 
Substitution into continuity gives a 2D elliptic equation for ^ . 

Substitution into the definition of streamwise vorticity gives a 2D 
elliptic equation for This is solved coupled with the transport equation 

for which replaces both transverse momentum equations. Anderson [114] 
uses the same splitting as Briley ana McDonald and solves for ^ , vv', 
and However, the Poisson equation for pressure b is added to the 

system. ' 

Approximation by Algebraic System 

All of the methods considered here use finite difference techniques to 
approximate the viscous equations with an algebraic system. First order 
upwind differences are used in the main flow direction and secund order 
central differences in the cross-plane. Within this common framework, a few 
variations deserve comment. Patankar and Spalding [108] use a staggered grid 
similar to that of Harlow and Welch [115]. Velocity components and pressures 
are stored at different locations within a grid cell in oroer to simplify the 
differencing of the convective terms. The other methods store all variables 
at common locations within the gria. For large transverse velocities both 
Patankar and Spalding [108] and Ghia et. al. [110] switch to upwind 
differencing in the cross-plane. Roberts and Forester [111] add explicit 
local damping to deal with this problem. The other methods have not as yet 
encountered this problem. Anderson [113] in his 2D method applied Keller's 
box scheme to a system of first order partial differential equations. In the 
current version [114] for 3D flows, he goes back to a system of second order 
equations and uses differencings similar to the other methods. 

Solution of Algebraic System 

Finally, we consider the techniques for solving the algebraic system of 
equations at each cross-section. Patankar and Spalding [108] use a completely 
non-iterative scheme. Provisional velocity components are obtained from^the 
momentum equations by an ADI marching technique. The bulk correction to is 
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obtained b> satisfyin9 a global mass balance over the cross-section. The 
local correction is obtained from the 2D elliptic equation b3> several ADI 
sweeps. All velocities are next corrected using the approximate relationships 
between velocity and pressure corrections. Finally, the energy equation is 
solved again by ADI. The procedures of^riley [109] and 6hia et. al. [110] 
are similar to the above, except that is determined iteratively with the 
streamwise velocity and the 2D elliptic equations for <^> and are solved by 
point SOR. The overall process still is essentially non-iterative at each 
cross-section. Roberts and Forester [111] follow a sequence similar to that 
of Briley, without the equation for . However, the entire sequence is 
repeateo iteratively at each cross-section using updated pressures in the 
momentum equations until convergence is achieved. Briley and McDonald [112] 
solve for the streamwise velocity, density, and total enthalpy with a coupled 
block ADI technique and iterate to determine . Next scalar ADI is used to 
find the potential ^ . Finally, a coupled block ADI scheme is used to 
find Uf and while satisfying the no-slip condition at the wall. Anderson 

[114] solves a fully coupled system for primary velocity, ()) , p , 

and total enthalpy. At present point SOR is used with a 6 x 6 block inversion 
at each point. 

Application of Methods 

The methods of this section have been applied to a variety of flows. The 
method of Patankar and Spalding p08] has been used for developing laminar 
flow in a square duct with a moving wall [108] and a round turning duct 

[116] . Briley [109] calculated the developing laminar flow in rectangular 
ducts and included the effects of transverse buoyancy. The method of Ghia et. 
al. [110] has been applied to the developing laminar flow in straight ducts of 
polar cross-section [110] ano to turning ducts of rectangular cross-section 

[117] . Roberts and Forester [111] computed the turbulent flow in a 
rectangular-to-rounci diffusing transition duct. The method of Briley and 
McDonalo [112] has been used for laminar flows in a turning duct similar to a 
turbine blade passage [112] and for turbulent flow in a rectangular turning 
duct [118]. Finally, Anderson [113] has computed several 2D turbulent flows 
in axisymmetric ducts with curved walls and Anderson and Hankins [114] have 
applied their 3D method to the hot turbulent flow in a turbofan forced mixer 
nozzle. 

To illustrate the capability of these methods we present results obtained by 
Kreskovsky, Briley, and McDonald [118] for turbulent flow in a rectangular 
duct with a 90° bena. Figures (35) and (36) show computed primary and 
racial velocity profiles, for a cross-section 77.5° around the bend, 
compared with the LDV measurements of Taylor et. al. [119]. Note the very 
large radial velocity near the suction side of the channel. The quantitative 
disagreement with the data may be due in part to the large streamwise step 
used in the computation. 

Two general comments would seem in order at this point. First, all of the 
methods in this section seem promising provided that small enough streamwise 
steps are used in the computation and sufficient pains are taken to ensure 
that local continuity is satisfied accurately enough for correct secondary 
velocities to form. Second, from the standpoint of efficiency a non-iterative 
procedure would seem to be preferable. 

Finally, we note that Baker and Orzechowski [120] have developed a finite 
element parabolic method. Like the methods of Briley and Ghia et. al., it 
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sol ves 2D elliptic equations for 4^ and ^ . Like that of Anderson it uses 
the small transverse velocity assumption to parabolize the system. As of this 
writing, however, we are not sufficiently familiar with the operation of the 
method to discuss it further. 


PARTIALLY PARABOLIC METHODS 
Satisfaction of Local Continuity 

Three of the four methods considered in this section are similar to fully 
parabolic methods and hence our discussion can be somewhat abbreviated. We 
first consider the satisfaction of local continuity, both Pratap and Spalding 
[121] and Moore and Moore [122] employ approximate relations between the 
velocity and pressure corrections and adjust these variables at each 
cross-section in the same manner as Patankar and Spalding. However, in 
successive passes, as the pressure field is refined, these approximate 
corrections should approach zero. Chilukuri and Pletcher [123] correct the 
velocity field at each cross-section through a potenti'al in a manner 
similar to that of Briley and Ghia et. al. By assuming to be zero at the 
downstream station the method simultaneously corrects the primary velocity 
without the need for a bulk pressure correction. As the pressure field is 
refined in successive passes 4^ should approach zero. Dodge [124] splits 
the velocity into viscous and potential parts by setting 

V = U 

In this expression U is obtainea by marching the momentum equations. (j)is 
updatea after each full sweep by substituting the above expression into the 
continuity equation and solving the resulting three-dimensional elliptic 
equation. Unlike the other methods, however, ^ does not approach zero with 
successive passes through the grid. 

Elliptic Pressure Update 

The technique for updating the elliptic pressure field is the main 
distinguishing feature of thes>i methods. Pratap and Spalding [121] use the 
pressure field obtainea from the continuity corrections during the march. An 
ad hoc means of distributing these corrections upstream is mentioned in the 
paper but not discussed. Moore and Moore [122] obtained an elliptic pressure 
correction equation from an approximation to the divergence of the vector 
momentum equation. This is solved after each march of the viscous equations 
using source terms calculated and stored during the march. The source terms 
and therefore the corrections approach zero after many marching passes through 
the field. Chilukuri and Pletcher [123] use the pressure Poisson equation 
obtained from the divergence of the full momentum equation. This is solved 
after each march with source terms again evaluated and stored during the 
march. The new pressure replaces the old. Dodge [124, 125] introduces an 
approximate inviscid relation between the pressure and the velocity 
potential ^ . Once the 3D elliptic equation is solved for 4^ , after each 
march, the pressure |> is obtained from the algebraic relation. 

Approximation by Algebraic System 

As in the previous section, finite difference techniques are used in each 
method to approximate the differential system by a large algebraic system. In 
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the marching equations first or second order upwind differences are used in 
the primary flow direction and second order central differences in the cross- 
plane. ine main distinguishing features are discussed briefly. Both Pratap 
and Spalding [121] and Chilukuri and Fletcher [123] use the same staggered 
grid employed by Patankar and Spalding. The velocity components and the 
pressure are each stored at different locations in a grid cell. Moore and 
Poore [122J use a different staggered scheme with all velocity components 
stored at common locations and only the pressure corrections stored at 
different places in the grid. Dodge [125] stores all variables at common grid 
locations. In addition he introduces finer subgrids near the walls in order 
to resolve the viscous layers. The elliptic pressure update equations of both 
Moore and Moore [122] and Chilukuri and Pletcher [123] use central differences 
in all directions. The global potential equation of Dodge [125] uses the 
mixed upwind-central differencing of reference [18]. 

Solution of Algebraic System 

The solution techniques for the algebraic systems are also similar to those of 
the preceding section. All methods use an ADI technique to obtain the 
velocity components from the momentum equations during the march. Since 
Chilukuri and Pletcher [123] solve 20 problems, they only need to perform a 
tridiagonal matrix inversion in one direction. Pratap and Spalding [121] and 
Moore and Moore [122] use ADI to solve the pressure correction equations at 
each cross-plane. Only the Moore's, however, iterate with the momentum 
equations until the pressure corrections are acceptably small. Both Moore and 
Moore [122j and Chilukuri and Pletcher [123] use point relaxation procedures 
to solve the elliptic pressure equations. The Moore's omit points in the near 
wall region of the boundary layer in their method. Dodge [124, 125] obtains 
separate marching solutions on each of his near wall subgrids and couples 
these to the interior marching solution at their common boundaries. The 
global potential equation for is solved by the transonic relaxation 
technique of reference [18], with the near wall subgrid points omitted from 
the field. Thus, Dodge is the only method able to compute viscous transonic 
flows with shocks. 

Application of Methods 

The partially parabolic methods have been applied to a somewhat broader range 
of flows than those of the previous section. The method of Pratap and 
Spalding [121] has been used for 3D turbulent flow in rectangular turning 
ducts [126]. Moore and Moore have computed 3D turbulent flow in an 
accelerating rectangular elbow [127] and two centrifugal impellers [128, 

129]. Chilukuri and Pletcher [123] computed 2D laminar flow in the inlet of a 
straight channel over a broad range of Reynolds numbers. Dodge [125] has 
calculated 3D turbulent flow in a rectangular diffuser and a low aspect ratio 
turbine stator. We show results from two of these computations as examples of 
the state-to-the-art. 

Stanitz et al . [130] measured the turbulent flow in an accelerating 
rectangular elbow designed by means of potential flow theory. The planform is 
shown in figure (37). Moore and Moore [127] computed this flow for cases with 
exit Mach numbers of 0.26 and 0.4 in which spoilers were used to thicken the 
incoming endwall boundary layers. Figure (38) shows computed and measured 
wall static pressure on four potential surfaces from inlet to exit. Figure 
(39) siiows computed and measured total pressure loss contours at the exit 
plane downstream of the bend. 
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Dodge [125] has computed 3D turbulent flow through a highly- loaded 
low-solidity high-aspect-ratio turbine stator. His results are compared to 
the measurements of Waterman [131]. Figure (40) sliows computed contours of 
static pressure ratio on the blade suction surface. Figures Hi) and (42) 
show computed and measured static pressures for the hub and tip blade 
sections, respectively. 

In both of these cases there is good qualitative and fair quantitative 
agreement between the computations and the measurements. These results are 
encouraging, but much work remains to be done. In particular, the development 
of an accurate and efficient global pressure update procedure for transonic 
viscous flows with shock waves would be a major accomplishment. 


ELLIPTIC METHODS 

All of the methods in this category are capable of computing separated flows. 
Some of them employ parabolic approximations in the viscous terms only. The 
utility of these methods in the analysis of separated flows justifies their 
inclusion in the present section. The methods which solve the compressible 
equations in conservation-law form have shock-capturing capability. This is 
not generally available in the parabolic methods to date. In the absence of 
separation or shock waves, the elliptic methods may or may not be more 
accurate than the parabolic ones. This is primarily dependent on the accuracy 
of the differencing of the inviscid terms in the equations. Most of the 
methods now in use for internal flows are adaptions of techniques discussed 
here. However, other methods, especially those now under development for 
external aerodynamic applications, will undoubtedly be adapted for internal 
flows in the near future. 

Methods for Steady Equations 

We consider first, techniques for solving the steady viscous equations. A 
popular method introduced by Caretto et al. [132] is based on that of Patankar 
and Spalding [109]. The method uses the so-called SIMPLE algorithm which 
stands for Semi Implicit Method for Pressure Linked Equations. In this 
technique an initial pressure field is substituted into the momentum equations 
which are in turn solved for provisional values of the velocities. An 
approximate velocity-pressure correction relation is then substituted into 
continuity to give a 3D equation for pressure corrections. The corrected 
pressure field, underrelaxed for stability, is then substituted into the 
momentum equations to continue the process. Many variants of this method have 
been studied by Raithby and Schneider [133]. They found that reintroduction 
of the time derivatives, with backward time differencing, into the momentum 
equations increased the convergence rate of the algorithm. Of course, this 
converts the method to an implicit time marching procedure. 

A different procedure has been used by Walitt et al. [134]. The 3D steady 
equations are transformed to a 2D unsteady system by treating one direction of 
the problem as time-like and evaluating the spatial derivatives in this 
direction from a previous solution. The equations are solved by marching 
through the field in the time-like direction. After one complete sweep, the 
time-like direction is switched and the equations are marched in a new 
direction. Since an explicit procedure is used to solve the equations on each 
cross-plane, very small steps in the marching direction are required in order 
to assure stability. This method has been applied to a centrifugal impeller 
[134] and to flow in a supersonic compressor cascade with splitter vanes [135]. 
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Methods for Unsteady Equations 

We next consider techniques for solving the unsteady viscous equations. These 
are closely related to explicit and implicit time marching methods for the 
Euler equations, although in some cases the viscous method came first. In an 
explicit technique all spatial derivatives are evaluated at a previous time. 
Hence, the solution proceoure is unchanged by the addition of viscous terms. 
The boundary conditions on the velocity are changed in order to enforce no- 
slip at solid walls. The close mesh spacing, necessary to resolve the 
boundary layers, imposes severe time step restrictions in order to maintain 
numerical stability (CFL limit). Bosman and Highton [136] have developed an 
explicit viscous method closely related to the Euler method of reference 
[79]. The intended application here was to 3D subsonic flow in rotating 
machinery. Shang et al. [137] have implemented a 3D version of MacCormack's 
explicit predictor-corrector scheme [71] on a vector computer. The 
application here was to supersonic shocked flow in a rectangular wind tunnel. 
Spradley et al. [38] have also implemented a 3D explicit method on a vector 
computer. The time updating was quite similar to the MacCormack scheme. 
However, the spatial discretization was obtained by me^ans of the general 
interpolants method (GIM), [139]. The application here was to supersonic flow 
in an exhaust nozzle. 

In an implicit technique, the non-linear spatial derivative terms are 
linearized about the previous time, and backward differencing is used on the 
time derivative. The resulting coupled linear system is modified by the 
addition of the viscous terms. However, for central differencing of both 
inviscid and viscous terms, the matrix structure remains the same as for the 
Euler equations. These implicit techniques permit high resolution of the 
viscous layers without severe time-step restrictions. Briley and McDonald 
[92, 93] and Beam and Warming [94] have developed similar implicit techniques 
which were described earlier. Briley and McDonald have concentrated on 
viscous flows from the outset. Steger [96] has implemented the Beam and 
Warming procedure in a 2D curvilinear-coordinate, viscous code, which is 
applicable to a broad range of flow conditions including those in 
turbomachinery. In addition, Steger drops the streamwise viscous diffusion 
terms in the equations. 

Ghia et al. [140] have developed what is termed a semi -el 1 iptic implicit 
method for 2D incompressible flow. They first present a fully elliptic method 
that solves the complete momentum equations together with a Poisson equation 
for the pressure. Continuity is enforced by driving one of the source terms 
in the Poisson equation to zero in the manner of Harlow and Welch [115]. The 
semi-elliptic method is obtained by dropping the streamwise viscous diffusion 
terms. This in turn permits a simpler solution procedure for the momentum 
equations. 

Approximation by Algebraic System 

All but one of the elliptic methods use finite difference techniques to 
approximate the differential equations by a large algebraic system. Spradley 
et al. l138] use the Glh formulation which is similar to the finite volume 
approaches discussed earlier. In this approach the dependent variables are 
represented by interpolating functions over the interior of local mesh 
volumes. The algebraic system is obtained from weighted integrals of the 
differential equation over each mesh volume. Only Caretto et al . [132] use 
the staggered grid differencing scheme discussed earlier in reference to 
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Patankar and Spalding [108]. Bosman and Highton [136] employ the staggered 
arrangement discussed earlier as part of their Euler method [79]. All the 
other methods store all of the dependent variables at common grid locations. 
The hybrid central-upwind differencing used by Patankar and Spalding has also 
been used by Caretto et al. [132] and Briley and McDonald [93]. Ghia et al. 
[140] has used upwind differencing of the streamwise convective term 
everywhere in the flow field. 

Solution of Algebraic System 

The techniques for solving the algebraic system of equations have been 
discussed earlier, either in connection with the associated Euler methods or 
in the two sections on parabolic methods. Only a few comments are made here. 
All of the explicit time marching procedures as well as the explicit spatial 
marching method of Walitt et al. [134] use local update schemes somewhat akin 
to point Jacobi relaxation. Values at the new time depend only on a few 
surrounding values at the previous time. This is a slowly converging 
process, but it is easily coded. All of the implicit time marching procedures 
as well as the steady method of Caretto et al. [132] use ADI techniques to 
solve at least the momentum equations. These techniques can be made very 
fast, especially if the time-step is permitted to vary both in time and space, 
see e.g. McDonald and Briley [141]. Finally, for completeness, we note that 
the elliptic pressure equation has been solved by ADI, Caretto et al. [132], 
and by point SOR, Ghia et al. [140]. 

Application of Elliptic Methods 

To illustrate the state-of-the-art for elliptic methods in internal viscous 
flows, we present results from four of the analyses discussed above. The 
actual range of applications of the methods is too broad to be covered here. 

Humphrey et al. [142] have computed laminar flow in a square turning duct 
using a method based on that of Caretto et al. [132]. Results were compared 
with LDV measurements by the same authors. Figure (43) compares calculated 
and measured (circles) streamwise velocity profiles at mid-span and 
quarter-span for several cross-sections progressing from upstream to 
completion of the 90° bend. Tl.e disagreement at the last three stations may 
be due to inadequate grid resolution (10 x 15) over the cross-section. A 
similar computation and comparison with data for turbulent flow in the same 
duct is presented in reference [143]. 

Buggeln et al. [144] nave used the method of Briley and McDonald [93] to 
compute both laminar and turbulent flow in curved ducts, channels and pipes. 

We present comparisons with the data of Taylor et al. [119] for the same 
turbulent flow case computed by Kreskovsky et al. [118] with the parabolic 
method of Briley and McDonald [112]. This also is the same duct used for the 
studies reported in references [142] and [143]. Figure (44) shows comparisons 
for streamwise velocity profiles at the symmetry plane for several stations 
upstream of, around, and downstream of the bend. Figure (45) shows 
comparisons for several radial velocity profiles at the 77.5° station. 
Comparing Figure (45) with Figure (36) shows considerable agreement between 
the two methods. 

Shang et al. [137] have computed the 3D flow in a square wind-tunnel diffuser 
for a case with a normal shock wave system interacting with the turbulent 
tunnel-wall boundary layer. Figure (46). Figure (47) shows computed Mach 
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number contours in the plane of symmetry. Figure (48) shows comparisons with 
the LDV measurements of Abbiss et al. [145] for the streamwise velocity in the 
symmetry plane through the interaction region. The agreement is quite good. 

Finally, Ghia et al. [140] have computed 2D laminar flow in a channel with a 
large constriction using both their fully elliptic and semi -el 1 iptic methods. 
Figure (49) shows computed streamline contours for a case with a large 
separation zone on the downstream side of the constriction. The upper part of 
the figure shows an enlargement of the central portion of the full channel 
shown below. Figures (50) and (51) show computed results of both methods for 
wall shear and wall static pressure, respectively. The agreement between both 
methods is excellent even for the wall shear in the separated region on the 
lower wal 1 . 

In general these methods all perform quite well. The disagreement with 
experiment where it exists is probably due to inadequate grid resolution. 
However, much work remains to be done to determine the numerical accuarcy of 
all these methods. For 3D flows this can be a very expensive process. 


TURBULENCE MODELING 

The recently completed AFOSR-HTTM-Stanford Conference on Complex Turbulent 
Flows was organized for the purpose of providing an assessment of the 
state-of-the-art in turbulent flow prediction, especially that of turbulence 
modeling. A wide range of test flows with reliable experimental data was 
assembled and computational groups were invited to submit computed results for 
these flows to the Conference. The comparisons of these results with the data 
together with the findings of the Evaluation Committee will be presented in 
reference [7]. We confine ourselves to a few observations on these 
proceedings. 

The current focus in turbulent modeling seems to have shifted from the one and 
two equation eddy viscosity models of a few years ago to methods which predict 
the Reynolds stresses themselves. Many results obtained with these models 
were presented at the Conference. However, while the promise of these models 
is considerable, at the present they show little or no advantage over the 
simpler treatments. This finding is especially true for separated flows. 

The results of any turbulent flow prediction depend at least as much on the 
numerical technique as they do on the turbulence model. For most of the flows 
included in the conference, it was not possible to separate the limitations of 
the numerics from those of the turbulence models. Grid refinement studies 
were presented in only a few cases and many of these were inadequate. Hence, 
in judging the merit of any turbulent flow computation, the method had to be 
considered as an amalgam of numerical procedure and turbulence model. 

In many internal flows, and especially in those through turbomachinery 
components, as a result of the rapid turning of the fluid, the evolution of 
the flow is dominated by the balance between strong pressure gradients and 
centrifugal forces. Turbulent mixing, though present, plays a lesser role. 
Indeed the development of complex secondary flows in these components, while 
dependent on the presence of shear layers due to upstream viscous effects, is 
an essentially inviscid phenomenon. Hence, if one can develop a sufficiently 
accurate and efficient numerical technique for solving the partial differential 
equations of fluid dynamics, good predictions of such flows should be possible 
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even with a relatively simple turbulence model. This view is supported by the 
observations of Humphrey et al. [143] in their studies of turbulent flow in 
turning ducts. 


CONCLUDING REMARKS 

The current status of potential flow methods for turbomachinery is quite 
advanced. A large number of codes have been developed for analyzing 
two-dimensional transonic flow with shocks on blade-to-blade surfaces. While 
an order of magnitude speed up in computing time may be possible with the 
latest procedures, many of these codes run fast enough for routine use by 
designers. In addition, we are aware of three codes which can compute 
three-dimensional potential flow in rotors. One of these has demonstrated 
transonic shock-capturing ability. Most of the stream function codes for 
turbomachinery applications are based on classical relaxation or matrix 
inversion techniques. They solve for subcritical flow on blade-to-blade or 
hub-to-shroud surfaces. Recently, however, a transonic stream function method 
has been developed and one application to a cascade is under development. In 
the case of the primitive variable Euler equations a large number of codes 
have been developed for turbomachinery applications. Most of these solve for 
two-dimensional, blade-to-blade flow with shocks. There are, in addition, at 
least three codes which analyze three-dimensional shocked flows in rotors. 

With few exceptions, however, these codes utilize older long-running numerical 
algorithms whose accuracy is not up to current standards. Since solution of 
the Euler equations is now one of the most active areas of research in 
computational fluid mechanics and since a number of promising new methods are 
under development, this situation should be much improved in a few years. 

The state-of-the-art for viscous methods is much less developed than for 
inviscid ones. A few single-pass parabolic marching codes have been developed 
for three-dimensional flows in ducts and turning passages. These methods are 
relatively fast and fairly accurate in the absence of strong secondary flows. 
When strong secondary flows are present at high Reynolds numbers, either 
excessive grid or implicit numerical dissipation may be needed to maintain 
stability in inviscid regions of the flow. A few multi-pass partially 
parabolic codes have been developed for three-dimensional flows in turning 
passages and at least two of these have been applied to turbomachinery 
rotors. In these methods the pressure field is treated as elliptic and is 
updated as the computation proceeds. Much work remains to be done to speed up 
these methods and to extend the pressure correction techniques into the 
transonic regime. A number of codes have been developed which solve the full 
time-averaged Navier-Stokes equations in either two or three dimensions, in 
both ducts and turbomachinery, for both subsonic and transonic flow. The 
comments on accuracy and computational efficiency which were made for the 
Euler equation methods are equally valid here except that the computer times 
are even longer. It is likely that the status of Navier-Stokes codes will be 
improved substantially with the advent of new improved Euler solvers since 
most regions of viscous flows are dominated by inviscid effects. Finally, we 
come to the question of turbulence models for all the viscous methods. Here 
we can close on a more hopeful note. Although it is not true that turbulence 
models are getting much more accurate, in many flows in turning passages and 
turbomachinery the velocity field is determined primarily by the balance 
between centrifugal forces and pressure gradients and the effects of 
turbulence are relatively weak. Hence, if one has a sufficiently accurate 
numerical procedure, one can hope to adequately compute such flows even with a 
relatively simple turbulence model. 
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FIG. 7 - 0-TYPE GRID FOR COMPRESSOR BUDE ROW 
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F16. 11 - DULIKRAVICH - NACA 0012 COMPARISON WITH CAUGHEY 
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fig. 13 - EFFECT OF STREAf'TUBE THICKNESS AND RADIUS CHANGE ON FLOW 
ABOUT SUPERCRITICAL STATOR HUB SECTION 



FIG, lA - HIRSCH - VKI TURBINE CASCADE COMPARISON WITH EXPERIMENT 
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FIG. 17 - BCSIWN - LOW SPtED COMPRESSOR WITH MERIDIONAL GRID 



FIG. 18 - BOSNAH - COMPRESSOR BLADE SHAPE WITH SI MEAN STREAMLINES 
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FIG. 19 - RELATIVE VELOCITY DISTRIBUTIONS ON HUB WITH NON-ITERATED 
AND ITERATED BOSTAN SOLUTIONS 



F16. 20 - RELATIVE VELOCITY DISTRIBUTIONS ON SHROUD WITH NON-ITERATED 
AND ITERATED BOSMAN SOLUTIONS 



FIG. 21 - HIRSCH - IMPELLER WITH FINITE ELEMENT GRID AND STREAMLINES 
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- HIRSCH - CALCULATED AND EXPERIMENTAL STATIC PRESSURE 
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FIG. 25 - HIRSCH - CALCULATED MACH NUMBER CONTOURS ON IMPELLER SHROUD 



FIG. 26 - MEASURED MACH NUMBER CONTOURS NEAR TIP 
OF NASA TRANSONIC COMPRESSOR ROTOR 





FIG. 27 - CALCULATED NACH NUflBER CONTOURS NEAR TIP OF NASA TRANSONIC 
COMPRESSOR ROTOR - THOHPKINS METHOD 
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FIG. 28 - BOSMAN - RADIAL INFLOW TURBINE HUB WITH CALCULATED PATHLINES 





FIG. 29 - BOSMAN - RADIAL INFLOW TURBINE BLADE SURFACE': 
WITH CALCULATED PATHLINES 
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FI6. 31 - BROCHET - PRESSURE CONTOURS IN PLANE OF SW.ETRY 
FOR CONVERGING SIDEWALL COMPRESSOR CASCADE 
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FIG. 55 - BRILEY - COMPUTED AND MEASURED AXIAL VELOCITY PROFILES 
AFTER 77.5° TURNING IN CIRCULAR ARC SQUARE DUCT 
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FIS. 56 - BRILEY - COMPUTED AND MEASURED RADIAL FLOW VELOCITY PROFILES 
AFTER 77.5° TURNING IN CIRCULAR ARC SQUARE DUCT 
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FIG. 37 - STAN1T2 RECTANGULAR ELBOW WITH STREAMLINES AND POTENTIAL LINES 
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FIG. 38 - MOORE - CALCULATED AND MEASURED WALL STATIC PRESSURE ON SUCTION, 
PRESSURE AND END HALL SURFACES OF STANITZ ELBOW 
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FIG. AO - DODGE - COMPUTED STATIC PRESSURE ON SUCTION SURFACE 
Of HIGHLY LOADED TURBINE STATOR 
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FIG. - DODGE - CALCULATED AND flEASURED STATIC PRESSURES 
FOR TURBINE STATOR HUB SECTION 
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FIG. A2 - DODGE - CALCULATED AND MEASURED STATIC PRESSURES 
FOR TURBINE STATOR TIP SECTION 
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FIG. « - HUMPHREY - CALCULATED AND MEASURED PROFILES AT MID-SPAN 
AND QUARTER-SPAN AROUND CURVED RECTANGULAR DUCT 
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FIG. W - BRILEY - COMPUTED AND MEASURED AXIAL VELOCITY PROFILES 
ALONG SYMMETRY PLANE OF CURVED SQUARE DUCT 
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FIG. US - BRILEY - COrPUTEB AND MEASURED RADIAL FLOW VELOCITY PROFILES 
AFTER 77.5° TURNING IN CIRCULAR ARC SQUARE DUCT 



FIG. A6 - SUPERSONIC WIND TUNNEL WITH NORMAL SHOCK WAVE 




FIG. 48 - SHAfiG - COMPUTE) AND ftASUREI VELOCITY IN SYWETRY PLANE 
THROUGH INTERACTION REGION 
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FIG. 49 - GHIA - COfiPUTED STREAMLINE CONTOURS FOR CHANNEL KITH CONSTRICTION 
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FIG. 50 - GHIA - COMPARISON OF SEMI-ELLIPTIC AND ELLIPTIC COMPUTATIONS 
FOR HALL SHEAR 
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